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1.  Introduction 

Practical  engineering  calculations  of  flow  fields  in 
which  boundary  layer  separation  occurs  have  historically 
been  very  difficult  with  realistic  methods  appearing  only 
recently  [1  -7].  The  main  difficulty  in  making  such  calcu- 
lations lies  in  the  fact  that,  even  for  high  Reynolds 
numbers,  the  classical  boundary  layer  theory  is  not  appli- 
cable. Although  the  effects  of  fluid  viscosity  and  rotation 
may  indeed  be  confined  to  a (limited)  region  of  the  flow, 
typically  such  "boundary  layers"  are  thick  and  induce  sig- 
nificant deflection  of  the  flow  streamlines  in  the  vicinity 
of  the  boundary  surfaces.  This,  of  course,  causes  the 
pressure  distribution  on  solid  boundaries  to  be  different 
from  the  "inviscid"  pressure  distribution  and  hence  the 
pressure  must  be  treated  as  an  unknown  in  both  the  main  body 
of  the  flow  and  in  the  "boundary  layer"  region.  Additional 
complications  are  introduced  into  the  problem  by  the  possible 
presence  of  turbulence  and  compressibility  effects. 

Tnere  are  three  possible  approaches  one  might  adopt  in 
attempting  to  predict  flows  exhibiting  boundary  layer  sepa- 
ration. These  are: 

1.  Solve  the  (time  averaged)  Navi er-Stokes  equations 
including,  if  necessary,  turbulence  models  of 
varying  sophistication. 

2.  Use  a form  of  vi scous -i nv i sc i d interaction  theory 


in  which  the  flow  is  divided  into  (at  least)  two 


2 


regions,  in  one  of  which  viscous  forces  are  neg- 
lected. The  pressure  is  treated  as  an  unknown  in 
the  "boundary  layer"  eq^Jations  and  is  determined 
via  an  interaction  between  the  viscous  and  inviscid 
regions. 

3.  Introduce  "engineering  approximations"  in  which  the 
details  of  the  separated  flow  are  not  computed  but 
"forced"  on  the  flow,  e.q.  by  specifying  a fixed 
dividing  streamline  shape  and  treating  it  as  a solid 
boundary . 

Each  method,  of  course,  has  its  advantages  and  drawbacks. 
Method  1 has  been  the  subject  of  much  research  but  solutions 
obtained  to  date  have  tended  to  concentrate  mainly  on  shock- 
induced  boundary  layer  separation  or  separation  induced  by 
abrupt  changes  in  geometry  such  as  cavities  or  steps. 

Method  3 has  seen  much  use  in  engineering  practice  but 
the  necessary  experimental  information  always  restricts 
application  of  this  method  to  those  situations  which  closely 
duplicate  the  one  on  which  the  model  is  based. 

Method  2 adopts  the  middle  ground  in  that,  although 
lacking  complete  generality,  it  should  have  the  ability  to 
reveal  significant  details  of  the  flow.  In  method  2 (as  well 
as  3)  research  can  be  concentrated  on  basic  "modules"  (e.g. 
the  viscous  layer,  the  inviscid  flow,  the  interaction  mechan- 
ism) independently.  Previously  obtained  knowledge  and  exper- 
ience in  boundary  layer  calculations,  for  example,  can  be 
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brought  to  bear  because  the  viscous  layer  equations  are 
still  parabolic.  Within  the  framework  of  the  strong  inter- 
action approach,  many  levels  of  sophistication  are  possible. 

It  is  of  utmost  importance  that  the  problem  be  formulated 
correctly,  with  the  mechanism  of  interaction  clearly  de- 
lineated. Once  this  has  been  achieved,  the  analysis 
methods  in  each  region  can  be  chosen  with  some  freedom. 

For  example,  it  is  possible  to  use  either  integral,  finite 
difference,  or  finite  element  formulations  for  the  boundary 
layer  equations,  while  using  completely  different  methods 
for  the  inviscid  flow  region. 

In  the  research  described  in  this  reoort,  an  attempt 
was  made  to  develoo  a method  for  calculating  flows  in  which 
turbulent  boundary  layer  separation  occurs.  The  flow  fields 
were  restricted  to  those  in  which  the  maximum  Mach  number 
is  less  than  1.  Although  the  methods  investigated  have  many 
portions  that  are  applicable  to  both  plane  2-dimensional 
and  axisymmetric  geometries,  major  emphasis  is  placed  on 
axisymmetric  geometries.  The  strong  interaction  method  was 
adopted  as  the  basis  for  the  work.  Separate  development  of 
methods  for  calculating  separating  turbulent  boundary  layers 
with  necessary  “free  stream"  information  input  given  and  for 
calculating  inviscid  flows  with  appropriate  strong  interaction 
type  boundary  conditions  was  undertaken. 

Although  a complete  calculation  was  not  obtained,  sig- 
nificant success  in  the  development  of  each  component  was 
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achieved.  The  failure  of  the  overall  method  is  believed  to 
be  due  to  the  lack  of  the  ability  to  make  an  initial  "guess'* 
close  enough  to  the  final  answer.  (Each  component  was 
developed  and  tested  with  "exact"  experimental  information 
input  as  the  necessary  information  from  the  other  component.) 

It  is  felt  that  the  use  of  an  interactive  computer  system 
(presently  unavailable  to  the  authors)  would  open  the  possibil- 
ity of  making  complete  calculations. 

In  the  remainder  of  this  report,  the  problem  will  be 
formulated  and  the  approaches  considered  in  this  research 
will  be  described.  First  the  problem  formulation  will  be 
discussed,  with  major  emphasis  on  the  interactive  mechanism 
adopted.  Then  the  calculations  for  each  flow  region  will 
be  described.  Finally  the  attempts  at  integration  of  a 
complete  model  will  be  described  and  suggestions  for  future 
research  made. 


n 

2.  Problem  Formulation 

Consider  the  separated  flow  configuration  shown  in  Fig- 
ure 2.1.  For  some  portion,  the  boundary  layer  is  attached 
to  the  body.  Because  of  adverse  pressure  gradients,  the  flow 
separates  from  the  body  and  a region  of  reversed  flow  (sepa- 
ration bubble)  develops.  At  some  downstream  location,  the 
separation  bubble  is  terminated  by  either  a reattachment  to 
a solid  surface  or  a realignment  of  the  flow  via  interaction 
with  a separated  boundary  layer  originating  on  the  "lower" 
side  of  the  body.  It  is  here  assumed  that  this  separated 
region  is  relatively  "thin"  and  that  the  flow  is  steady  i.e. 
vortex  shedding  and  bluff  body  flows  are  excluded.  This 
type  of  separation  often  occurs  near  the  trailing  edge  of 
airfoils  and  off  of  boattailed  after  bodies  employed  on  jet 
engi  ne  nacelles  and  missiles.  In  the  latter  case,  the  reattach- 
ment may  in  fact  be  onto  a propulsive  plume,  although  in  the 
work  considered  herein,  reattachment  onto  a solid  surface  was 
assumed.  It  is  assumed  that  the  Mach  number  of  the  flow  is 
everywhere  less  than  1 and  that  the  Reynolds  number  is  suf- 
ficiently large  that  boundary  layers  are  turbulent.  All 
surfaces  are  assumed  to  be  smooth  and  impermeable. 

In  analyzing  this  flow,  the  major  assumption  that  is 
made  is  that  the  flow  can  be  divided  into  two  distinct  regions: 

(1)  a rotational  "viscous"  region  in  which  all  of  the  effects 
of  shear  stress,  turbulence,  reverse  flow,  etc.  are  concen- 
trated and  (2)  an  inviscid,  irrotational  outer  region.  The 
viscous  region  is  assumed  to  be  finite  in  extent  and  within 
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FIGURE  2.1  SCHtXATIC  OF  SEPARATED  FLOW 
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this  region,  it  is  assumed  that  the  standard  boundary  layer 
assumptions  regarding  the  neglect  of  stream-wise  stress 
gradients  is  applicable.  Thus  the  governing  equations  for 
the  viscous  layer  are 

Conti nui ty 


IM  + l£V  + i 
3x  3y 


pu  dR 
R 17 
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Momentum 


(2.1) 


3x  3y  " p dx  p 3y 


(2.2) 


In  equation  2.2.  it  is  assumed  that  the  pressure  is 
constant  across  the  viscous  layer;  this  assumption  may  be 
subsequently  replaced  with  a "centrifugal  correction"  of 
the  form: 


or  in  integral  methods  of  solution  with  an  assured  poly- 
nominal  P(y)  distribution  with  a typical  parameter  (say 
P(6)  - P(o))  determined  from  solution  of  the  (integral)  y- 
momentum  equation  [8  - 9 ]. 

It  Is  here  emphasized  that  it  is  not  necessarily  assumed 

that  P(x)  is  known  a priori  for  use  in  equation  2.2;  at  this 

point  it  represents  one  of  the  unknowns.  If  P is  assumed 

uniform  across  the  viscous  layer,  then  we  may  replace  it 

with  P (x).  The  thickness  of  the  viscous  region  is  denoted 
e 

by  6. 
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When  dealing  with  compressible  flows,  the  energy 
equation  must  also  be  considered.  In  this  work,  all  compres- 
sible flows  were  assumed  to  be  adiabatic  and  a solution  to 
the  energy  equation  was  provided  by  employing  the  Crocco 
integral  relating  the  temperature  and  velocity  in  the  viscous 
1 ayer . 


• t + RF  ( ^ ) M|  (I  - ♦M 


(2.3) 


The  recovery  factor  was  computed  from 
RF  = Pr'^^ 


Outside  of  the  viscous  layer,  the  flow  is  assumed  to  be 
irrotational  and  is  governed  by 

3(PV)  * j . 0 (2.4) 
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The  pressure  in  the  outer  region  is  governed  by  the 
Euler  equations 
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particular  at  the  boundary  between  the  two  regions 


e ^ 


du 

^"e  dT 


(2.6) 


9 


The  boundary  conditions  which  can  be  applied  directly 


u « 0 ^ 
V * 0 J 


at  solid  boundaries 


(2.7) 
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The  interaction  between  the  two  regions  is  introduced 
into  the  problem  via  the  matchi ng  condi tions  which  must  be 
satisfied  at  the  boundary  between  the  two  regions.  If 
y = 5(x)  is  the  thickness  of  the  viscous  layer,  then  the 
matching  condition  is  that  the  velocity  (vector)  at  y = o(x) 
must  be  the  same  for  the  two  regions.  This  is  most  easily 
accomplished  by  equating  the  total  velocity  (or  Mach  number) 
and  the  flow  angle  at  the  edge  of  the  viscous  layer  to  the 
corresponding  values  from  the  inviscid  region.  Thus  we 


wri  te 


(v5  + vM^/2  » (ui  + vi) 


Ml  > M I , 

e'outer  e'inner 


9 = O + a 


(2.9) 


5(x) 


(2.10) 


Equations  2.9  and  2.10  provide  the  formal  connection  between 
the  two  regions. 

It  should  be  noted  here  that  vi scous  - i nvi sci d inter- 


action is  usually  accounted  for  via  "displacement  thickness 


interaction"  rather  than  a formal  matching  of  the  flow 
fields  at  the  edge  of  the  viscous  layer.  In  displacement 
thickness  interaction,  it  is  assumed  that  the  velocity  (or 
angle)  at  the  edge  of  the  viscous  layer  can  be  obtained  by 
calculating  the  flow  over  the  boundary  surfaces  augmented 
by  adding  the  boundary  layer  displacement  thickness.  This 
is  an  approximation  to  equations  2.9  and  2.10  in  that 
the  angle  O is  very  closely  the  angle  of  the  displacement 
surface.  Displacement  thickness  interaction  has  been  used 
in  most  other  ( non-supersonic ) separated  flow  analyses 
[1-5];  however,  in  this  study,  the  more  formal  interaction 
of  equations  2.9  and  2.10  has  been  retained. 

With  the  introduction  of  the  customary  equations  of 
state  for  ideal  gases,  equations  2.1  - 2.10  provide  a closed  set 
sufficient  to  determine  the  flow.*  It  is  not  possible  to 
solve  the  equations  simultaneously,  however,  so  an  inter- 
itive  technique  is  employed,  with  the  viscous  and  inviscid 
regions  computed  alternately.  The  matching  conditions  then 
appear  as  boundary  conditions  for  one  region  or  the  other, 
with  the  most  current  information  from  the  alternate  calcu- 
lation being  used. 

In  the  sections  that  follow,  methods  for  making  calcu- 
lations in  each  region  will  be  discussed.  In  developing 
the  analysis  for  each  region,  it  is  assumed  that  complete 
information  on  the  conditions  at  the  inner-region  boundary 
is  available  i.e.  either  or  both  of  M^,  ©are  known.  First 

* For  compressible  flow,  the  adiabatic  energy  equation  is 
included  in  the  inviscid  region. 
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3.  Viscous  Reflion  Calculations 

In  developing  a calculation  method  for  the  viscous 
region  ("boundary  layer")  with  separation,  it  is  assumed 
that  the  necessary  information  from  the  inviscid  region 
(pressure,  velocity,  or  Mach  number  and/or  flow  angle  at 
the  edge  of  the  boundary  layer)  is  available.  Only  one  of 
these  variables  is  needed  as  input  and  in  fact  in  a full 
interaction  calculation,  only  one  can  be  provided;  the  other 
one  being  calculated  by  the  analysis  itself.  In  the  attached 
boundary  layer  (weak  interaction)  region  it  is  the  pressure 
that  is  input,  with  the  flow  angle  being  calculated.  At 
this  point  it  is  necessary  to  assume  that  if  either  of 
these  variables  is  needed,  it  can  be  provided.  The  task  of 
the  inviscid  region  analysis  then  will  be  to  compute  the 
necessary  variable  which  is  not  calculated  in  the  boundary 
layer  analysis. 

The  task  then  becomes  to  construct  a solution  to  equa- 
tions 2.1  and  2.2,  given  the  conditions  of  equation  2.7  and 
one  of  the  conditions  2.9,  2.10.  Historically,  three  methods 
have  been  developed  to  solve  these  equations.  These  are  (1) 
integral  methods,  (2)  finite  difference  methods,  and  (3)  finite 
element  methods.  Although  each  of  these  methods  has  several 
points  in  its  favor,  as  well  as  several  drawbacks,  it  was  de- 
cided at  the  outset  of  this  research  to  develop  an  integral 
method,  based  on  the  following  reasons: 


w 
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(1)  Integral  methods  were  felt  to  be  simpler  to 
program, 

(2)  It  was  felt  that  an  integral  method  should  consume 
less  computer  time, 

(3)  An  integral  method  can  be  formulated  which  solves 
directly  for  the  parameters  necessary  to  formu- 
late the  viscous  - inviscid  interaction  i.e.  6(x), 

Mg,  0,  while  in  a finite-difference  or  finite 
element  method,  parameters  such  as  5(x)  and  hence  0 
are  ill-defined, 

(4)  By  selecting  a velocity  profile  which  has  an 
inherent  reverse  flow  capability,  regions  of 
reverse  flow  can  be  calculated  while  still  "marching" 
in  the  downstream  direction,  without  the  necessity 

of  employing  a special  algorithm  for  reverse  flow 
points , 

(5)  Local  inaccuracies  in  turbulence  models  may  be 
smoothed. 

As  a result  of  this  study,  it  is  possible  to  draw  some 
conclusions  regarding  these  points;  this  will  be  dealt  with 
in  a later  section. 

3.1  Derivation  of  Integral  Equations 

In  employing  the  integral  method,  the  boundary 

layer  equations  2.1  and  2.2  are  multiplied  by  u'’’y”  and 

formally  Integrated  across  the  layer  from  y » 0 to  y = 5 
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(actually,  rather  than  integration  of  2.2  as  is,  2.1 
is  multiplied  by  u"’^^y^  and  added  to  2.2,  then  the 
result  is  integrated).  This  produces  a set  of  (usualfy 
independent)  equations  which  can  be  solved  for  a finite 
number  of  parameters.  Using  m=n  = 0 produces  the 
momentum  integral  equation,  m=l,  n*0  produces  the 
mechanical  energy  equation,  while  m»0,  n»l  produces 
the  moment  of  momentum  equation.  The  formal  applica- 
tion of  this  procedure  leads  to  the  following  equations 
following  lengthy  algebraic  manipulations 


Continuity  Equation: 


5 Ug  6 ‘•^0  ax  ^0  Ug  3x 


1-5-^- 

ex  e e 


(3.1  ) 


Momentum  Equation: 


u^  3 


I-  (^)d?  . 2 ^ I-  (MdS] 


0 2 3x  0 

U Mg 

e 


a OgUg  3x  'Ug 


^ ^ ^ 1 fiul.  j dC-1]  ^ ^ 

R dx  0 2 e “ 0 u'  '"e 

e e e e 


11  tan  0 
25  " 5 


(3.2) 


HHK'a 


Moment  of  Momentum  Equation; 


(H"  It  (|-)5d5  ♦ 2 / 


1 ou  3 / u 


0 u . 3x  ' 0 
e ^e 


fr  (^)5dC 


" Pjd,  3x 


0- ^ -4  ff]  ^cd5 

' '’e“e 


du 


<['-'1^  S-dT^^isf)  ^^§V^d]d5 


e e 


tan  Q 1 1_  .1  T 

5 2u.  6 0 2 

e Pe% 


dC 


(3.3) 


Mechanical  Energy  Equation; 


e Ug  -e 


du. 


1 pu  ’ 


[({3-M^}  ^ ^ + i P)  — dC 

e u„  dx  R dx ' 0 3 ^ 

e e 


_ du 

- 2_  e u_  d^j  , . tan_0 
Ug  dx  ' 0 Ug  p 


1 2t  3 /U  \ j,. 

6 *^0  r H 

e 

e e 


(3.4) 
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In  the  above  equations,  isentropic  flow  external 
to  the  boundary  layer  has  been  assumed  and  the  approxi- 
mation of  equation  2.6  has  been  employed. 

In  order  to  make  any  calculations  with  these  equa- 
tions, it  is  necessary  to  evaluate  the  integrals  appear- 
ing in  them.  This  in  turn  requires  knowledge  of  the 
variation  of  the  velocity  ratio,  density  ratio,  and  their 
X - deri vati ves  across  the  boundary  layer,  as  well  as 
the  variation  of  the  shear  stress  across  the  boundary 
layer.  The  usual  approach  is  to  assume  a velocity 
profile  of  the  form 

(x)  “ ^(y»a  j (x)  ,a  Jx) . . . ) 

where  the  parameters  a^  are  functions  of  x only.  The 
density  ratio  can  be  related  to  temperature  ratio  and 
hence  to  velocity  ratio  via  eq.  2.3.  The  shear  stress 
distribution  is  related  to  velocity  distribution  via 
Newton's  law  of  viscosity  and/or  a "turbulence  model". 

In  some  (if  not  most)  methods  in  actual  use,  the  shear 
stress  integrals  themselves  are  postulated  as  functions 
of  parameters  of  x.  Using  the  assumed  profile  the 
indicated  integrations  on  ? can  be  performed  and 
the  partial  derivatives  become  total  derivatives  with 
respect  to  x of  the  a^(x).  The  a^(x)  thus  effectively 
become  the  unknowns  of  the  problem  and  as  many  equations 
are  required  as  there  are  a^  to  calculate.  The  boundary 
conditions  2.7,  2.9,  2.10  are  usually  satisfied  auto- 


i 

\ 

matically  by  the  function  chosen  for  the  velocity  j 

profile;  the  problem  then  becomes  an  initial  value  prob-  | 

i 

1 em  for  the  a . ' s . 

1 

3.2  Choice  of  Velocity  Profile 

It  is  here  noted  that  henceforth,  consideration  is 
restricted  to  turbulent  flows.  When  chosing  the  velocity 

1 

profile,  the  following  restrictions  and  guidelines  must 
be  met 

(1)  The  profile  must  accurately  represent  the 
actual  velocity  distribution, 

(2)  The  profile  must  be  capable  of  exhibiting  re- 
verse flow  near  the  wall  (negative  wall  shear), 

(3)  The  parameters  involved  in  the  profile  should 
be  meaningful  and  convenient. 

Relevant  to  the  third  point,  it  is  noted  that  the 
integral  equations  themselves  (3.1  -3.4)  contain  the 
parameters  Ug(x)  (equivalent  to  fig(x)  for  isentropic 
flow  external  to  the  boundary  layer),  6(x),  C^(x),  and 
Q(x).  Now  Ug(x)  and  0(x)  are  exactly  what  is  needed  for 
the  interaction  with  the  external  flow.  In  addition, 

5(x)  must  be  known  so  that  the  boundary  between  the  two 
regions  may  be  located.  It  therefore  seems  that  Ug(x), 

5(x)  and  C^(x)  are  logical  choices  for  the  parameters 
of  the  velocity  profile.  Now  for  incompressible  turbu- 
lent non-separated  boundary  layers,  it  is  well  known  that 

L - , ^ 


the  wall-wake  velocity  profile  [10]  provides  an  excel- 
lent approximation  to  the  actual  velocity  profile. 

Alber  and  Coats  [11  ],  Mathews  et  al.  [12],  and  Alber 
etal.  [13]  have  shown  that  this  profile  can  be  modi- 
fied to  account  for  compressibility,  and  Alber  et  al  . [13] 
and  Kuhn  and  Nielsen  [2-3]  have  shown  that  the  profile 
can  be  modified  to  include  both  reverse  flow  and  a 
laminar  sublayer.  Accordingly,  the  following  form  is 
adopted  for  the  velocity  profile 


4>  = ^ j sin  {a\[i  In  (l+y'")  + B - (l.Sy'^  + B)e'-''®^^ 


+ aug  sin^ ( J ^) } 


(3.5) 


where 


k = .41  , B = 5.0 


a = 


1 + R(X|i)M^ 


. y"  =|A| 


C^/2 


C./2I 


1/2 


The  parameters  in  this  equation  are  u^ (x ) [M^ ( x ) ] , 

5(x) , A(x) , Ug(x) . 

The  unit  in  the  In  and  the  linear-exponential  term 
provide  a smooth  decrease  to  zero  velocity  at  y * y^  » 0. 

The  constants  in  the  linear-exponential  term  were  select- 
ed to  provide  an  optimum  fit  with  the  Spal di ng  - Kl ei ndi enst 
law  of  the  wall,  as  oresented  by  White  [15]. 


With  constant  pressure  across  the  boundary  layer 


assumed,  the  density  ratio  can  be  calculated  by 

^ [1 +RF 1 -(})^ ) ] (3.6) 

If  the  flow  is  incompressible,  p/p  =1  , v *v 

e we 

and  the  velocity  ratio,  4),  is  given  by  the  term  enclosed 
in  the  {}  in  Eqn.  3.5,  with  the  "a"  removed. 

Although  there  are  4 parameters  in  (3.5),  they  are 
not  completely  independent  since  the  condition 

u/u^  =1  at  y = 5 
e 

must  be  satisfied.  This  results  in 


j sin  In  (1+6^)  +B 


(1.55'^  + 8)e''^®‘^  ^ ^ 


F(Mg,X,5,Ug)  * 1 


(3.7) 


This  may  be  called  a "skin  friction  law",  at  any 
rate  it  is  a relationship  between  X{x) , 6(x),  Mg(x), 

Ua(x)  which  must  be  satisfied. 

P 

3.3  Formulation  of  Differential  Equations  for  X,6, 

^e’  *^8’  ® 

Within  the  framework  of  the  integral  approach, 
complete  information  on  the  viscous  layer  is  represented 
by  knowledge  of  the  five  parameters  M^,  O,  X,  5,  Ug. 

Now  one  of  these  parameters  (either  or  O)  is  specified 
by  the  inviscid  flow  coupling  so  4 (differential)  equa- 
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One  equation  that  must  be  satisfied  is  3.7,  the  other  3 
necessary  equations  can  be  selected  from  3.1  - 3.4.  The 
choice  is  actually  between  the  mechanical  energy  equa- 
tion (3.4)  and  the  moment  of  momentum  equation  (3.3), 
with  the  continuity  and  momentum  equation  being  re- 
tained in  all  cases.  For  the  sake  of  generality,  both 
equations  will  be  further  developed  here  although  the 
moment  of  momentum  equation  was  ultimately  selected  for 
reasons  to  be  discussed  later. 

Substitution  of  the  assumed  velocity  profile  and 
the  resulting  density  profile  3.5,  3.6  into  3.1  - 3.4 

allows  (at  least  in  theory)  the  performing  of  the  integrals 
and  the  differentiation  with  respect  to  x.  The  form  of 
the  resulting  equations  is 


§.  + A ^ + A + A §.  = A + A + 

1 dx  ''2  dx  ^3  dx  dx  ^5  ^6 


Where  the  coefficients  involve  M^,  X,  5,  Ug 

and  require  evaluation  of  an  integral  on  y from  0-^6 , A 

6 

involves  ©and  A^  is  possibly  a shear  stress  integral. 

Now  the  complicated  form  of  3.5  for  <<)  makes  analytical 
integration  of  and  its  powers  and  derivatives  as  needed 
in  3.1  -3.4  impossible  (at  least  for  compressible  flow), 
therefore  the  integration  must  in  general  be  performed 
numerically  to  find  A^-^A^.  Accordingly,  the  equations 
3.1  - 3.4  are  here  presented  with  the  integrals  yet  to  be 
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evaluated.  It  is  emphasized  that  although  the  integrands 
are  algebraically  complicated,  they  are  functions  of 
X,  6,  Ug  and  C only  and  hence  the  indicated 
integration  on  ^ can  be  performed  (numerically  if  neces- 
sary) if  M^,  X,  6,  Ug  are  known.  The  equations  are 

Continuity ; 
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IX  Tr^2  iT: 


(continued) 
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Moment  of  Momentum: 
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1 /■  _L_ 

6 ■*  0 2 

e e 


Mechanical  Eriergy  Equation 


(3.10) 


* fr  tS-H*’  If  * 3r*"]  * ♦>  |g-)d5 

e 3y  e e 

T , dM 

*[(3-M^)/;r^>d5  - 2/>d;]  dir 

0 e 

* If  * 3r*ndC)  ff 

3y 

* H [♦’  If  * ^ 

* IJr  [♦■  !?*  3r*nd5)  ^ 


tan  O 2 ,1  T 3<j)  j..  j dR 

“ - ^ 6 ^0  H ■ R d7  'I? 

(3.11) 

The  various  partial  derivatives  of  the  <|)  and  r(=*  p/pg) 
are  tabulated  in  Appendix  A. 

Equations  3.8,  3.9,  either  of  3.10  or  3.11  and  3.7, 
together  with  a specification  of  either  Mg(x)  or  ©(x) 
and  a set  of  initial  values  thus  form  a closed  set. 

Now  3.8-3.11  are  di  fferential  equations , which  will 
ultimately  be  solved  numerically,  while  3.7  is  an  algebraic 
equation.  In  performing  numerical  computations,  it  has 
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been  found  easier  to  solve  4 differential  equations 
than  3 differential  and  one  algebraic  equation;  there- 
fore, 3.7  is  differentiated  with  respect  to  x to  yield: 


it  _ rlF  ^ 3F_  ^ + fit  + it_  iili  ^ 

dx  '■3a  3M^  ..+  3M  ^ dx  '•3X  3X  ^ dx 

e do  e do 


. r 3F_  361t  ^ ^ - n 

‘•gg+  35  ■>  dx  3Ug  dx  ■ 


(3.12) 


The  partial  derivatives  appearing  in  the  above  are  tabu- 
lated in  Appendix  A. 

We  have  thus  completely  formulated  the  viscous 
layet"  problem,  with  the  exception  of  a means  of  evalu- 
ating the  shear  stress  integrals  appearing  in  the 
mechanical  energy  and  moment  of  momentum  equations. 

This  of  course  requires  the  introduction  of  some  semi- 
empirical  "model"  of  the  turbulent  transport  process; 
this  complicated  subject  will  now  be  dealt  with. 


3.4  Turbulence  Models 


In  order  to  close  the  set  of  equations  3.8,  3.9, 
3.10  or  3.11,  and  3.12,  it  is  necessary  to  evaluate 
either 


/!  -4-  dC 
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for  the  moment  of  momentum  equation,  or 
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for  the  mechanical  energy  equation.  For  turbulent  flow, 
T represents  the  sum  of  the  laminar  viscous  shear  and 
the  Reynolds  shear  stresses. 


3u  — 
U gy  - puv 


■^L 


Complete  information  on  the  Reynolds  stress  is 
not  currently  available,  hence  these  terms  (or  their 
integrals)  are  evaluated  via  some  semi -emperi cal  model. 

In  many  integral  methods,  an  attempt  is  made  to  intro- 
duce information  on  the  integrals  themselves  [ 1 4 ] ; 
however,  in  differential  methods  a local  formulation 
for  T is  needed.  In  the  current  research,  since  a 
numerical  evaluation  of  integrals  of  the  velocity  pro- 
file function  and  its  derivatives  (typical  terms  appear- 
ing on  the  right  hand  sides  of  equations  3.8-3.11)  is 
to  be  undertaken,  a numerical  integration  of  a local 
shear  stress  formulation  was  selected.  An  eddy  viscosity 
formulation  was  selected  in  order  that  a negative  shear 
would  be  predicted  in  regions  of  reverse  flow. 

T - (u  + pc)  (3.13) 

It  is  well  known  that,  turbulence  is  a phenomenon 
which  exhibits  at  least  two  length  scales,  therefore  for 
turbulent  flow  near  walls,  separate  formulations  for 
the  eddy  viscosity  near  the  wall  and  in  the  outer  region 
are  often  used.  A popular  choice  for  the  eddy  viscosity 
is  the  mixing  length  formulation 
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e » (3.14) 

The  mixing  length  is  at  least  an  order  of  magni- 
tude larger  in  the  outer  region  of  the  boundary  layer 
than  near  the  wall.  Most  eddy  viscosity  models  employ 
a mixing  length  formulation  near  the  wall;  however,  in 
the  outer  region  an  alternate  form  may  sometimes  be 
used,  typical  is  the  Clauser  model 


^outer 


.016  u^5. 
e k 


(3.15) 


Both  approaches  have  been  applied  to  a wide  variety 
of  turbulent  boundary  layer  calculations,  including 
separating  flows  [2  - 4];  in  the  basic  ("equilibrium") 
form  embodied  in  3.14  and  3.15  there  seems  to  be  little 
to  choose  between  the  two.  Calculations  undertaken  in 
the  current  research  bear  this  out  for  attached  boundary 
layers;  however,  the  mixing  length  model  was  found  to  be 
more  ameable  to  extension  to  non-equilibrium  flows 
and  to  boundary  layers  with  backflow;  therefore,  the 
mixing  length  formulation  was  selected  as  the  basic 
turbulence  model  for  the  present  work.  The  mixing  length 
typically  exhibits  the  following  characteristics  [is] 


Mear  the  wall:  l-y 

3/2 

Very  near  the  wall:  il-y  ' 

In  the  outer  region:  t-i ndependent  of  y,  t-S 


The  following  (continous)  mixing  length  distribu- 
tion embodies  all  of  these  characteristics  and  was  employed 


in  this  work 


I 

5 
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exp 


r.  _:_  (-W  + 1 IP  y i 


1/2 


2 6v 


p dx 


3) 


(3.16) 

The  last  term  represents  the  van  Driest  [16]  damp- 
ing factor  as  modified  for  the  effect  of  pressure  gradient 
as  suggested  by  Cebeci  et  al.  [17].  For  "equilibrium" 
attached  turbulent  boundary  layers,  the  value  of  is 

a constant  equal  to  aporoximately  0.09. 

Now  when  a boundary  layer  separates  from  a solid 
surface,  the  damping  effects  of  the  wall  on  the  turbu- 
lence are  removed;  in  fact,  the  layer  becomes  almost 
entirely  "wake  like",  accordingly  the  following  model  was 
adopted  for  the  mixing  length  in  separated  boundary  layers. 
With  reference  to  Figure  3.1,  the  layer  is  divided  into 
two  regions  at  the  zero  velocity  point  y^ . Above  the 
zero  velocity  point 


Z = = ( 2.^/5  ) 6 (3.17) 

while  below  the  zero  velocity  point 

i (3.19) 

Having  introduced  the  mixing  length  model  the 
shear  stress  integrals  may  now  be  formulated.  Assuming 
a laminar  viscosity  law  of  the  form 

M/Ug  • (T/Tg)^/"^  (3.20) 


the  shear  stress  function  becomes 


« (A)  liii 
u^6  '3c’ 


Substituting  the  mixing  length  equation  3.16  and 
introducing  the  velocity  profile  parameters,  we  obtain 


R M j 7/4 

iTT  " (r)  * 
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(a|x|-c|-5^)  ])  iff, 


(3.22) 


for  the  attached  boundary  layer  and 


2 » , 
.3^1 

^6  ' 'sc' 


— = < 
Ug5  ^ 


(3.23) 


The  shear  integrals  appearing  in  the  moment  of 
momentum  equation  (3.10)  and  mechanical  energy  equation 
(3.11)  can  be  evaluated  as  functions  of  the  velocity 
profile  parameters  M_,  X,  6,  u-  and  their  derivatives 

6 W 

using  3.21,  3.22,  3.23,  3.5,  and  the  derivatives  relations 
of  Appendix  A,  if  l/S  is  specified  (t  /6  becomes,  in 
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effect,  an  additional  parameter  of  the  flow  in  the 
boundary  layer). 

3.4.1  Non-Equilibrium  Effects  in  the  Turbulence  Model 

For  equilibrium  attached  boundary  layers,  the 

value  of  I /5  is  rather  well  determined  as 
00 

IJ6  = 0.09 

essentially  independent  of  pressure  qradient  effects. 
Successful  calculations  of  non-equilibrium  boundary 
layers  have  also  been  made  using  this  value;  however, 
the  simple  use  of  a constant  value  for  this  parameter 
was  found  to  be  unsatisfactory  for  calculations  of  the 
separating  boundary  layers  considered  here.  At  least 
two  reasons  for  this  behavior  may  be  offered: 

(1)  "History"  effects  are  present,  especially 

near  separation.  As  separation  is  approached, 
the  boundary  layer  grows  rapidly,  that  is  6 
increases  rapidly.  The  turbulence  structure 
does  not  respond  rapidly  so  that  as  5 grows 
rapidly,  does  not  and  thus  decreases, 

(2)  The  separating  boundary  layer  is  in  reality 
a process  of  transition  between  an  attached 
shear  flow  and  a free  shear  flow.  Even  for 
equilibrium  layers,  while  a value  of  l^/S  = .09 
is  appropriate  for  attached  flows,  a value  of 
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3 0.05  to  0.07  is  appropriate  for  free 
shear  layers  [18]-  A calculation  method 
for  both  attached  and  free  shear  flows  and 
the  transition  between  them  (separating 
boundary  layers)  must  reflect  this  change  in 

00 

In  order  to  successfully  calculate  separating 
boundary  layers  using  the  current  method,  it  is  neces- 
sary for  to  change  as  the  flow  develops  toward 

separation,  separates,  develops  a free  shear-like  charac- 
ter, and  (possibly)  reattaches  and  redevelops  as  an 
attached  boundary  layer.  It  is  of  course  necessary  to 

connect  the  change  of  2.^/5  with  changes  in  their  flow 

★ 

parameters.  There  are  two  possibilities  for  accomplish- 
ing this: 

(1)  Use  an  algebraic  formulation  of  the  form 

= f(Mg,A,6,Ug) 

(2)  Solve  a differential  equation  of  the  form 

iljs)  • »_/«.«) 

Within  the  differential  equation  approach,  the  choice 
is  between  solving  an  empirical  differential  e.quation  of 
the  "lag"  or  "departure  from  equilibrium"  type 


* The  idea  of  a formulation  I /&  - prescribed  function 

00 

of  X is  rejected  as  being  possible  only  if  the  answer 
is  known  in  the  first  place! 
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d(Z^) 

dx 

or  of  solving  a differential  equation  with  a more 
solid  physical  basis,  typically  some  form  of  the 
"turbulence  kinetic  energy  equation". 

During  the  course  of  this  investigation,  all 
three  of  the  above  approaches  were  investigated. 

Each  approach  will  be  described  in  turn. 

Considering  first  the  algebraic  formulation,  the 
overall  approach  of  Alber  [19]  was  employed.  The 
reasoning  proceeds  as  follows.  For  equilibrium 
boundary  layers,  Clauser  [20]  has  shown  that  the  ap- 
propriate correlating  parameter  for  turbulence  para- 
5* 

meters  is  S(=  ~ ^)  thus  for  equilibrium  boundary 
^w  ax 

layers  we  would  assume 
= f(S) 

The  function  f must  satisfy 

f = 0.09  6=0  (flat  plate) 

f -*•  0.05  to  0.07  as  6 ^ (free  shear  layer) 

The  form  chosen  here  was  determined  by  numerical 
experimentation  and  is  typical  of  other  methods  [ 2 ] ; 
it  is 

I 

* 0.055  + 0.035  exp  (-8/5)  (3.24) 

* In  terms  of  the  principal  parameters  employed  in  the 
present  study  * 2 ^e  1 ‘^^e 
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Following  Alber  [ 19]  it  is  assumed  that  the 
equilibrium  correlation  can  be  used  if  6 is  related 
to  the  equilibrium  shape  parameters  rather  than  its 

it 

definition  as  a pressure  gradient  parameter  . Accord- 
ingly, 6 is  calculated  from 

G = 6.1  /e  + r.Sl"  - 1.7  (3.25) 


where  G (the  equilibrium  shape  parameter)  is  related 
to  the  local  velocity  shape  parameters  by 
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(3.26) 


Equations  3.24  - 3.26  are  sufficient  to  relate 
£ /5  to  the  local  velocity  profile  parameters  M ,A, 

00  g 

d.Ug. 

Plotted  in  Figures  3.2  -3.4  are  skin  friction  coef- 
ficients, boundary  layer  thicknesses,  and  shape  factors 
predicted  for  the  incompressible  two  dimensional  sepa- 
rating boundary  layer  flow  of  Simpson  et  al . [ 21  ] . 

For  reasons  that  will  become  clear  later,  only  the  re- 
sults up  to  = 0.0005  are  shown.  Clearly,  the  predic- 
tions with  the  constant  (equilibrium)  value  of  £„/5  are 
unsatisfactory  (in  fact  it  is  not  even  possible  to 
predict  separation  using  this  model).  In  fact,  predic- 
tions with  the  algebraic  correlation  are  good  enough 
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In  Alber's  terminology,  3 is  "unhooked"  from  the 
pressure  gradient. 
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that  a differential  equation  method  would  have  to  be 
significantly  better  or  at  least  be  more  general  in 
order  to  justify  the  added  complexity  (and  computer 
time)  of  solving  an  additional  differential  equation. 

The  simplest  differential  equation  which  may  be 
employed  to  predict  turbulence  is  the  empirical  "lag" 
equation  apparently  first  suggested  by  Goldberg  [22] 
and  employed  successfully  in  the  turbulent  boundary 
layer  methods  of  Nash  and  Hicks  [23]  and  White  [15]. 
For  mixing  length,  this  equation  must  be  written 


^ ^ const  s 
dx  6 ^ “■’eq 


(3.27)  i 


] 


where  I ) = 0.09  5 . 

“ eq 

The  value  of  the  constant  in  the  equation  is  on 
the  order  of  0.1  and  can  be  determined  by  numerical 
experimentation.  It  is  important  to  note  that  the 
equation  cannot  be  written  in  the  form 


dx 


const 

6 


( 0.09 


because  if  at  any  point  (say  an  initial  point  far 

upstream  of  separation)  = 0.09,  it  will  always 

retain  this  value!  In  terms  of  I /6,  equation  3.27 

00 

becomes 


” /_!_  d6  X 

6 ^6  dx' 


[0.09 


(3.28) 
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Calculations  for  the  Simpson  et  al  . [21]  flow 
and  for  the  two  dimensional  transonic  separating  flow 
of  Alber  et  al  . [13]  were  preformed  using  3.28  and 
were  not  significantly  different  from  those  using 
3.24  - 3.26. 


The  differential  equation  that  has  been  used  most 

extensively  to  represent  turbulence  is  the  "turbulence 

kinetic  energy  equation".  This  equation  has  been 

employed  in  differential  [24  ] and  integral  [25,26] 

forms  and  has  been  used  with  differential  [24,26]  and 

integral  [25]  methods  of  boundary  layer  predictions. 

In  the  current  research,  the  basic  approach  of  McDonald 

et  al.  [25,26]  was  employed  to  convert  the  turbulence 

kinetic  energy  equation  into  a differential  equation  for 

I /6.  The  development  is  as  follows.  The  partial  dif- 
00 

f erenti a 1 equation  governing  the  turbulence  kinetic 
energy  is 
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(3.28) 


where  $ is  a collection  of  many  fluctuating  terms  and 
may  be  found  elsewhere  [ 15].  The  terms  in  the  equa- 
tion are  usually  given  names  of  the  following  form: 

Advection  ■ shear  stress  turbulence  production  + 
normal  stress  turbulence  production  - pressure  - 
strain  diffusion  - viscous  dissipation  (of 
turbul ence ) 
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Equation  3.28  may  be  written  as  an  integral  equa- 
tion by  adding  q*  x (continuity  equation)  and  inte- 
grating on  y from  0 -*■  6 . Assuming  negligible  free  stream 
turbulence,  the  result  is 


- dy  - J/f  ^ dy 


(3.29) 


Following  McDonald  et  al.  [25,26],  Bradshaw  et  al. 

[ 24  ],  and  Collins  and  Simpson  [27]  it  is  assumed  that 


■ Fiji 


§7 

r ^ 

^t  3y 


_ Normal  Stress  Production 
Shear  Stress  Production 


(3.30) 


T 3/2 

<t  = (Fiji)  /L 


where 


a = const  = 0.15 
1 

F » function  of  x only  a 


1 

' ■'  max 
shear 


L “ dissipation  length  » L(y/6) 

It  is  pointed  out  that  only  Collins  and  Simpson 
consider  normal  stress  turbulence  production  significant; 
the  models  of  McDonald  et  al.  [25,26  ] and  Bradshaw  et  al. 
[ 24  ] are  obtained  by  setting  F = 1.  Substituting 


3.30  into  3,29  and  expres'iing  the  resulting  equation  in 
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terms  of  the  parameters  M ,a,6,Uq,  we  get 

6 D 

0 n o A 


dx 
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where 


T = 


't  _ _t_  ^ 

» u ' ■ ''e* 

e e 


(3.32) 


If  the  mixing  length  model  (equations  3.22  and  3.23) 

is  used  and  if  F is  formulated  in  terms  of  the  assumed 

velocity  profile  function  and  its  derivatives,  then 

Equation  3.31  becomes  a differential  equation  involving 

fi«,  X,  5,  Uo,  I /6,  and  their  derivatives;  in  effect  a 
e jj  00 

differential  equation  for  (It  still  remains  to 
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specify  L/6  » f(C)).  It  must  be  noted  that  the  formu- 
lation is  still  not  complete,  since  all  of  the  partial 
derivatives  of  |t|  have  not  been  presented.  They  will 
not  be  written  here  but  the  following  two  points  are 
1 i sted 

(1)  The  partial  derivatives  involve  many  terms, 
due  to  the  algebraic  complexity  of  the  mixing 
length  model  itself, 

(2)  It  is  necessary  to  evaluate  all  of  the  second 
derivatives  of  the  velocity  profile  function. 

After  all  operations  have  been  performed  and  the 
indicated  integrations  performed  (numerically)  equations 
3.8,  3.9,  3.10  or  3.11,  3.12  and  3.31  are  five  coupled 
ordinary  differential  equations  for  the  six  parameters 
0,  M , X , 5 , u , 2.  /5  and  can  be  solved  if  either  of  <H? 

0 (30  00 

or  Mg  are  specified.  Calculations  were  performed  using 
this  set  of  equations  for  the  Simpson  et  al.  [21  ] and 
Alber  et  al.  [13  ] flows.  The  results  were  as  follows. 

It  was  first  assumed,  following  McDonald  et  al. 
[25,26  ] Bradshaw  [24],  and  particularly  Collins  and 
Simpson  [27]  that  the  dissipation  length  is  a function 
of  y/<S  only  and  as  approximated  by 

? - r * 5.5  5‘] 

oo' 

for  attached  boundary  layers,  (3.32) 


for  separated  boundary  layers 
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where  L /6  = 0.09. 

00 

Calculations  using  this  formulation  were  quite 
discouraging,  being  similar  to  those  made  with  = 0.09, 

as  indicated  by  the  skin  friction  plots  of  Figure  3.5.. 

It  was  concluded  that  the  dissipation  length  formula- 
tion of  3.32  was  not  sufficiently  general. 

The  measurements  of  Simpson  et  al.  [21]  indicate 
that  in  a separating  boundary  layer,  the  dissipation 
length  (as  well  as  the  mixing  length)  in  the  outer 
portion  of  the  boundary  layer  decreases  as  separation 
is  approached,  accordingly  calculations  with  L^/6  = 0.07 
were  made;  these  were  inaccurate  in  regions  upstream 
of  separation. 

Realizing  that  L /6  in  fact  varies  with  the  flow, 
calculations  were  then  made  with 

L„/6  = 0.05  + 0.04  exp  (-6/4) 

These  calculations  were  slightly  better  but  by  now  it 
should  be  obvious  that  specifying  L^/6  is  a function  of 
the  mean  flow  parameters  is  no  different  than  specifying 
4^/5  as  a function  of  the  same  mean  flow  parameter  and 
that  the  latter  is  considerably  more  efficient!  As  a 
result,  it  was  concluded  that,  within  the  framework  of 
the  present  method,  there  is  no  advantage  to  using 
either  the  extremely  complicated  (algebraically) 
turbulence  kinetic  energy  equation  or  the  (artificial) 

"lag"  equation  and  the  simple  algebraic  mixing  length. 


a 
*.J1 


correlation  embodied  in  Equations  3.24  - 3.26  was 
adopted  as  both  simpler  and  more  accurate. 

3.4.2  Effects  of  Turbulent  Normal  Stresses 

Since  the  work  of  Newman  [28],  it  has  been 
recognized  that  as  a boundary  layer  approaches  sepa- 
ration, the  effects  of  the  turbulent  normal  stresses 
(pu^,  pv^)  are  not  negligible.  These  stresses  con- 
tribute the  term 

+ 1^  (pv^-cu^) 

to  the  right  hand  side  of  the  partial  differential 
momentum  equation  of  the  boundary  layer  (eqn.  2.2). 

In  addition  to  their  effect  on  the  mean  flow  equations, 
Simpson  and  co-workers  [21,27]  found  that  the  effects 
of  normal  stresses  on  turbulence  production  were  signi- 
ficant. Collins  and  Simpson  [27]  subsequently  presented 
a turbulence  kinetic  energy  equation  based  turbulence 
model  which  accounted  for  the  normal  stress  production; 
this  is  essentially  the  model  of  equations  3.30  - 3.32. 

The  work  involving  the  turbulence  kinetic  energy 
(integral)  equation  described  herein  thus  took  normal 
stress  turbulence  production  into  account. 

Returning  to  normal  stresses  in  the  mean  flow 
equations,  it  is  necessary  to  model  the  above  term  in 
terms  of  the  other  flow  parameters.  Following  Townsend 
[29],  McDonald  et  al.  [25,26  ] and  Collins  and  Simpson, 
it  is  assumed  that 
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~ and  ~ 

thus  p(v^-uM  ~ pq^  - |t^| 
resulting  in 

p(v^-uM  =*  ^ I7  I (3.33) 

A value  of  a » 2 was  selected  in  accord  with  previous 

3 

workers  [25-27]. 

Upon  integrating  the  normal  stress  term  across  the 
boundary  layer,  the  term 

-s  h 

is  added  to  the  right-hand  side  of  the  momentum  equation 
3.2  (m=n“0),  moment  of  momentum  equation  3.3  (m=0, 
n»l),  and  mechanical  energy  equation  3.4  (m  = 1 , n«0). 
Substituting  the  assumed  velocity  profile  and  the  mixing 
length  turbulence  model  and  expanding  results  in  the 
addition  of  the  following  terms  to  the  working  Equations 
3.9-3.11: 

dri 

To  the  coefficient  of  is  added: 

dC  + (2-Mp  ^ /;p'"c"lTldO 

To  the  coefficient  of  ^ is  added: 

,1  .m^n  3 I t|  j, 

S 5 3X 


^ 


-i.m  n 3|t| 

^ dS 


Attention  is  directed  to  the  required  partial 
derivatives  of  the  turbulent  shear  stress  function  Ixj. 
As  in  the  turbulence  kinetic  equation,  these  are  alge- 
braic functions  of  considerable  complexity,  involving 


second  derivatives  of  the  velocity  profile  function  with 
respect  to  the  4 parameters  M^,  X,  5,  Ug. 

Using  the  present  method  calculations  were  performed 
in  which  the  effects  of  normal  stresses  in  both  the  mean 


flow  and  turbulence  kinetic  equations  were  included.  As 
regards  the  turbulence  kinetic  energy  equation,  the  effects 
of  normal  stresses  were  not  found  to  be  significant  (the 
effect  of  normal  stresses  are  represented  by  the  "F" 
defined  in  Equation  3.30,  with  F = 1 when  normal  stress 
effects  are  neglected)  and  do  not  influence  the  conclu- 
sion of  the  previous  section  that  use  of  the  turbulence 
kinetic  energy  equation  itself  is  not  warranted.  As 
regards  the  effect  of  normal  stresses  on  the  mean  flow 
equations,  calculations  of  C^,  5,  H for  the  entire 
Simpson  et  al.  [21]  flow,  both  with  and  without  the 
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normal  stress  terms  were  carried  out.  No  difference 
in  predicted  values  of  C^,  5,  H,  Q,  and  were  observed; 
hence,  it  was  concluded  that,  within  the  present  frame- 
work, normal  stress  terms  do  not  greatly  effect  the 
accuracy  of  the  method. 

3.5  Mathematical  Details  of  Solution  Procedure 

In  terms  of  a closed  system  of  equations,  the 

boundary  layer  calculation  problem  is  now  completely 

formulated.  Complete  knowledge  of  the  boundary  layer 

★ 

flow  is  given  in  terms  of  the  parameters  0,  M^,  X, 

6,  u-  as  functions  of  x.  Via  the  velocity  profile 
specification  of  Section  3.  and  the  turbulence  model 
discussed  in  3.4,  Equations  3.8-3.12  become,  ultimately, 
five  first  order,  non-linear,  ordinary  differential 
equations  with  M^,  X,  6,  u-,0as  dependent  variables 

6 D 

★ ★ 

and  with  x as  the  independent  variable  . Now  informa- 
tion from  the  inviscid  outer  flow  must  be  input  to  the 
boundary  layer  method,  thus  either  Mg(x)  or0(x)  is 
regarded  as  a known  function  (the  possibility  of  a 
mixed  specification  with  M^  known  for  certain  x's  and  0 
known  for  other  x's  is  not  precluded).  To  calculate  the 
remaining  4 parameters,  only  4 equations  are  needed, 

* In  incompressible  flow,  is  the  ratio  of  the  local 
edge  velocity  to  a reference  velocity. 


♦♦  Using  a differential  equation  turbulence  model  intro- 
duces an  additional  parameter,  2.^/6,  and  an  additional 
differential  equation. 
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accordingly  either  the  moment  of  momentum  equation 
(3.10)  or  the  mechanical  energy  equation  (3.11)  is 
dropped.  The  ultimate  selection  of  one  of  these  over 
the  other  will  be  discussed  shortly;  the  other  critical 
question  of  whether  or  © is  taken  as  specified  will 
be  the  subject  of  a later  section. 

The  algebraic  complexity  of  the  velocity  derivative 
integrands  and  the  turbulent  shear  integrals,  as  well  as 
the  non-linearity  of  the  resulting  differential  equa- 
tions make  a numerical  solution  the  only  possibility. 

The  differential  equations  to  be  solved  have  the  form 

A.  -T— ^ + B.  ^ + C.  ^ + D.  = E . 

1 dx  1 dx  1 dx  1 dx  i 

i = 1 ,4 

The  coefficients  A,  B,  C,  D,  E are  functions  of  X, 

6,  Ug  and  involve  several  integrations,  the  coefficient 
E also  contains  ©.  The  integrations  in  A,  B,  C,  D can 
be  preformed  analytically  if  the  flow  is  incompressible 
[ 30];  however,  the  shear  stress  *i  ntegral  s in  E cannot 
and  for  compressible  flow,  neither  can  those  in  A,  B,  C, 
0.  In  the  present  method,  all  integrals  were  evaluated 
numerically.  A large  number  of  integrals  had  to  be 
calculated  (approximately  40  if  the  turbulence  kinetic 
energy  equati-on  is  to  be  solved,  30  if  not),  requiring  a 
rapid,  efficient  integration  method.  Because  of  its 
higher  accuracy  for  fewer  points  (that  is  fewer  calcu- 
lations of  the  integrand),  Gaussian  quadrature  formul as 
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were  selected  for  single  integrals.  The  double  inte- 
grals required  for  the  moment  of  momentum  equation  were 
evaluated  by  using  a 10  point  Simpson's  rule  for  the 
outer  integral  with  a Gaussian  quadrature  for  the  inner 
integral.  All  integrals  were  evaluated  in  10  strips 
from  C * 0 to  C = 1 . Because  of  the  more  rapid  varia- 
tions near  the  wall,  6 point  quadratures  were  used  for 
5 < 0.1,  4 point  for  0.1  < C < 0-3  and  2 point  for  ^ > 0.3. 

Because  of  the  appearance  of  the  same  terms  in  many 
of  the  integrands,  a great  computing  time  saving  was 
realized  by  calculating  integrals  simultaneously  rather 
than  one  at  a time. 

The  complications  involved  in  calculating  double 
integrals  would  seem  to  indicate  that  the  mechanical 
energy  equation  (3.11)  should  be  preferred  to  the  moment 
of  momentum  equation  (3.10).  Initial  efforts  were  in 
fact  concentrated  on  a method  using  the  mechanical  energy 
equation.  Numerical  experimentation  showed  that  extremely 
small  integration  steps  were  necessary  to  obtain  accuracy 
in  the  dissipation  integral 
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The  integrand  in  this  term  is  proportional  to  (9(i)/3C) 
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and  the  rapid  variation  of  3(J)/3S  near  the  wall  requires 
extremely  small  steps  of  integration.  For  this  reason, 
the  moment  of  momentum  equation  was  found  to  be  super- 
ior to  the  mechanical  energy  equation  and  was  adopted 
for  the  remainder  of  the  study. 

Initially,  a fourth  order  Runge-Kutta  method  was 
selected  for  solving  the  differential  equations;  how- 
ever, this  was  found  to  be  very  slow.  A fourth  order 
modified  predictor  - corrector  scheme  was  then  employed 
which  has  the  ability  to  change  step  size  to  speed  up 
calculations  or  increase  accuracy.  As  expected,  the 
predi ctor  - corrector  method  required  about  half  the 
computing  time  of  the  Runge-Kutta  method. 

In  order  to  start  the  calculations,  initial  values 

of  M , X,  (5,  Uq  must  be  known  at  some  x location.  When 
e D 

testing  the  program  using  experimental  data,  values  of 
Mg,  X(C^),  5 are  usually  available  at  a station  upstream 
of  the  region  of  primary  interest  and  calculations  are 
started  there.  The  value  of  Ug  at  the  initial  station 
is  obtained  from  the  "local  friction  law",  equation  3.7 
When  making  a complete  calculation  starting  from  a leading 
edge  or  front  stagnation  point,  theoretically  a laminar 
boundary  layer  calculation  and  transition  analysis  would 
be  needed;  however,  in  the  present  work  it  was  assumed 
that  transition  occurred  very  close  to  the  leading  edge 
and  the  calculations  were  initialized  by  setting  5=0 
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at  a leading  edge  and  using  turbulent  flat  plate  boundary 
layer  correlations  to  estimate  6 and  X(C^)  at  the  first 
point  downstream  of  the  leading  edge. 

3.6  Weak  Interaction,  Strong  Interaction,  and  the 

Appearance  of  Singular  Points 

The  coupling  between  the  viscous  and  inviscid  flow 
regions  requires  that  either  or  O,  as  calculated  by 
the  inviscid  flow  model,  be  specified  as  known  input  to 
the  boundary  layer  model.  Of  course,  in  classical 
boundary  layer  theory,  the  edge  velocity  (M^)  is  taken 
as  the  known  input  parameter.  This  is  referred  to  as 
the  weak  interaction  model  because,  to  first  order  ac- 
curacy for  these  boundary  layers,  the  presence  of  the 
boundary  layer  does  not  affect  the  pressure  distribution. 
In  this  approach,  the  (small)  effects  of  the  boundary 
layer  on  the  pressure  distribution  are  accounted  for  by 
adding  a "displacement  thickness"  to  the  body,  effectively 
modifying  the  body  shape.  If  the  weak  interaction  ap- 
proach is  adopted,  the  continuity  equation  (3.8)  can  be 
subtracted  from  the  momentum  and  moment  of  momentum  equa- 
tions (3.9,  3.10)  thus  removing©  as  a variable.  After 
these  modified  equations,  together  with  3.12  are  solved 
for  X,  5,  and  u^,  © can  be  calculated  from  the  contin- 
uity equation  (and  5*  can  be  calculated  by  integrating 
the  velocity  profile).  That  this  approach  to  the  calcula- 
tion of  all  types  of  boundary  layer  flows,  both  attached 


and  separating,  as  desirable  is  due  not  only  to  the 
simplification  resulting  from  the  elimination  of  © from 
the  system  of  equations  but  also  to  the  fact  that  methods 
to  calculate  inviscid  flows  with  body  geometry  speci- 
fied are  easier  and  considerably  more  plentiful  than 
methods  to  calculate  geometry  (i.e.  0)  with  velocity 
(Mg)  specified.  Unfortunately,  considerable  research 
has  demonstrated  that  the  weak  interaction  formulation 
is  inadequate  for  the  calculation  of  separating  boundary 
layers.  [1-5,31-34]  Typically,  aoolication  of  the  weak 
interaction  formulation  yields  one  of  two  results: 

(1)  If  a pressure  distribution  predicted  from  a 
completely  inviscid  analysis  of  the  flow  over 
a body  is  input,  the  weak  interaction  method 
typically  breaks  down  as  separation  is  approached, 
with  numerical  divergence  occurring.  Typically 

<S , 5*,  H,  and  © rapidly  increase  just  before 
breakdown.  The  occurrence  of  separation  can 
be  predicted  by  such  a method  and  its  location 
roughly  estimated,  but  no  flow  details  in  the 
vicinity  of  separation  can  be  obtained  and  the 
effects  of  the  boundary  layer  on  the  pressure 
distribution  cannot  be  obtained. 

(2)  If  an  experimentally  determined  pressure  distri- 
bution is  inout,  calculations  usually  do  not 
predict  separation;  instead  the  skin  friction 


levels  off  at  some  constant  value  near  the 
measured  separation  location.  Calculations 
by  Cebeci  et  al  . [3  2]  for  several  Incompres- 
sible flows  and  by  Gerhart  and  Bober  [33] 
for  the  compressible  flow  of  Alber  et  al.  [13] 
substantiate  this  conclusion.  Application  of 
the  present  method  (In  the  weak  Interaction 
mode)  to  the  Incompressible  separating  flow  of 
Simpson  et  al  [21]  produces  similar  results  as 
shown  In  Figure  3.2.  Cebeci  et  al.  and  Gerhart 
and  Bober  have  Indicated  that  It  Is  possible  to 
conclude  that  separation  Is  In  fact  occurring 
and  to  estimate  Its  location  from  such  calcula- 
tions, but  obviously  they  are  of  no  value  for 
predicting  details  of  the  flow. 

The  first  behavior  Is  not  a cause  for  major  concern 
because  In  a separating  flow  the  Inviscid  pressure  distri- 
bution Is  significantly  modified  by  the  v1  scous-invl scl d 
Interaction  so  some  calculation  scheme  which  artificially 
smooths  the  pressure  distribution  for  the  first  few  cycles 
of  Iteration  may  possibly  be  devised.  The  second  be- 
havior Is  more  decisive  In  Its  Implications  since  It 
Implies  that  even  If  the  exact  pressure  distributions  were 
somehow  arrived  at,  the  boundary  layer  flow  could  still 
not  be  accurately  calculated! 


The  I'ssolution  of  this  dilemma  has  been  presented 
by  a number  of  '•esearchers  [1-5,34]  and  requires  that 
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calculations  be  made  in  the  so  called  strong  interaction 
mode  in  which  the  pressure  is  considered  as  an  unknown 
in  the  boundary  layer  calculation  with  some  other  para- 
meter being  specified  as  input  information.  If  the 
interaction  formulation  of  the  current  method  is  followed 
(see  equations  2.9,  2.10),  the  obvious  variable  to  specify 
as  input  is  the  velocity  angle  at  the  edge  of  the  boundary 
layer©.  The  flow  angle  has  been  used  extensively  for 
coupling  in  supersonic  separated  flow  calculations  [35-37  ] 
but  has  not  been  employed  as  extensively  in  incnmpres- 
si bl e/subsoni c methods.  In  low  speed  flows,  the  dis- 
placement thickness  has  been  used  [3-5,38],  as  has  the  skin 
friction  coefficient  [2,3].  The  reason  for  chosing  one 
formulation  over  another  is  not  so  much  that  one  is  more 
appropriate  from  a physical  viewpoint  but  that  some  formula- 
tions may  be  more  easily  coupled  to  inviscid  flow  methods 
than  others.  In  this  regard,  there  is  no  way  that  initial 
guesses  of  C^(x)  can  be  iteratively  updated  by  solving 
for  the  flow  external  to  the  boundary  layer;  however, 
guesses  of  5*(x)  orQ(x),  which  lead  to  a M^(x)  calculated 
from  the  boundary  layer  equations,  can  be  updated  by  apply- 
ing the  calculated  Mg(x)  as  a boundary  condition  in  an 
inviscid  flow  analysis  (the  "inverse"  or  "design"  problem 
of  inviscid  flow  theory).  As  to  the  choice  between  6*  or 
0 as  the  strong  interaction  variable,  it  is  here  pointed 
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out  that jformal ly,  it  is  the  slope  of  the  velocity 
vector  on  the  surface  that  is  computed  by  an  inverse 
inviscid  flow  method,  with  the  corresponding  "body  sur- 
face" i.e.  body  plus  displacement  thickness,  determined 
by  integrating  the  surface  slopes.  In  addition, 

O naturally  appears  as  a variable  in  the  integral  boundary 

* 

layer  method,  6*  would  have  to  be  introduced  artificially. 
For  these  reasons,  it  was  decided  to  employ  ©as  the 
interaction  parameter  in  the  strong  interaction  mode  in 
the  present  work.  Therefore,  in  the  strong  interaction 
mode,  ©(x)  is  presumed  known  and  the  continuity,  momentum, 
moment  of  momentum,  and  skin  friction  equations  3.8  - 
3.10,  3.12  are  regarded  as  4 equations  for  calculating 
the  4 variables  X,  5,  u.. 

Initially,  it  was  proposed  to  calculate  the  entire 
boundary  layer  flow,  attached,  separating,  fully  separated, 
reattaching,  and  redeveloping,  using  the  strong  inter- 
action formulation.  This  was  ultimately  rejected,  it  being 


* This  is  most  easily  accomplished  by  putting  the  boundary 
layer  integral  equations  in  weak  interaction  form  to 
el imi nate  O and  introducing  the  definition  5*  » 

dy  =•  6*(M^,X,fi,uJ  » 5*(x)  as  an  extra 

equation,  which  may  be  retained  in  its  algebraic  form 
or  converted  to  a differential  equation  [3  ] . 


found  that  the  strong  interaction  equations  are  ill 
conditioned  in  regions  where  the  boundary  layer  is  non- 
separated.  It  is  therefore  necessary  to  calculate  the 
boundary  layer  flow  using  a combination  of  weak  and 
strong  interaction  formulations.  With  reference  to 
Figure  3.6t  calculations  are  initially  begun  in  the 
weak  interaction  mode,  near  a leading  edge  or  at  some 
point  well  upstream  of  separation.  Calculations  are 
marched  downstream  in  the  weak  interaction  mode,  until 
at  some  point  it  becomes  necessary  to  switch  to  the  strong 
interaction  mode  as  separation  is  approached.  Calcula- 
tions are  carried  through  separation  to  fully  separated 
flow  in  the  strong  interaction  mode.  In  flows  exhibiting 
large  separated  regions,  in  which  the  flow  essentially 
develops  as  a free  shear  layer,  it  may  be  necessary  to 
switch  back  to  the  weak  interaction  mode  to  make  calcula- 
tions in  the  free  shear  layer.  If  and  when  the  separated 
boundary  layer  reattaches,  it  once  again  becomes  neces- 
sary to  employ  the  strong  interaction  formulation  to  cal- 
culate through  the  reattachment.  At  some  point  downstream 
of  reattachment,  the  weak  Interaction  formulation  may  again 
be  used. 

The  obvious  question  that  arises  is  when  to  switch 
from  one  mode  of  calculation  to  the  other.  In  practice 
there  are  two  criteria  used: 
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(1)  When  separation  or  reattachment  is  approached 
calculations  must  be  switched  to  the  strong 
interaction  mode;  when  reattachment  is  "complete" 
calculations  must  be  switched  to  the  weak  inter- 
action mode. 

(2)  When  singular  points  (the  set  of  equations  being 
used  becomes  indeterminant)  are  approached,  it 
is  necessary  to  switch  to  the  opposite  formu- 
lation to  continue  calculations. 

From  previous  discussion  of  the  behavior  of  the  weak 
interaction  equations,  it  is  obvious  that  this  formula- 
tion must  be  abandoned  before  (fully  developed)  separation. 
The  nearest  to  separation  that  these  equations  might  pos- 
sibly be  extended  is  the  point  of  intermittent  separation, 
discussed  by  Sandborn  and  Kline  [3  9].  This  criterion 
indicates  that  intermittent  separation  begins  at  the  point 

where  , 

6 * 

> 1 Ml  - jJi) 

This  is  typically  satisfied  if  is  approximately  2.4. 

It  was  found  that  better  accuracy  could  be  obtained  by 
switching  even  before  this  criteria  is  satisfied;  therefore, 
the  current  method  used  two  criteria  to  indicate  the 
approach  to  separation  and  the  need  to  switch'  from  the 
weak  interaction  to  the  strong  interaction.  These  criteria 
are 
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> 2.0 

or  (3.34) 

O > .5 


The  former  criteria  tends  to  dominate  if  experi- 
mentally determined  inputs  (f1g(x))  as  specified  and  hence 
presumably  in  the  later  stages  of  iteration,  while  the 
latter  tends  to  dominate  if  ideal  flow  pressure  distri- 
butors are  used  (i.e.  in  the  early  stages  of  iteration). 
In  flows  redeveloping  following  reattachment,  calcula- 
tions are  switched  from  the  strong  back  to  the  weak 
interaction  mode  if 
"k  i 2-2 

and  (3.35) 


The  occurance  of  singularities  in  the  boundary 
layer  equation  is  a well  known  fact  the  stagnation 
point  and  separation  singularities  being  two  well  known 
examples.  Integral  boundary  layer  calculation  methods 
seem  to  be  especially  prone  to  exhibit  singularities, 
especially  in  the  vicinity  of  separation  and  reattachment. 
Singularities  may  be  divided  into  two  types,  those  which 
are  connected  with  significant  physical  occurrences  in 
the  flow  itself  and  those  which  are  not  necessarily  con- 
nected with  physical  occurrences  but  are  rather  attribu- 
table to  mathematical  anomalies  in  the  particular  calcu- 
lation method  chosen.  An  example  of  the  former  type  of 
singularity  is  the  well  known  Crocco-Lees  [31,35-37] 
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critical  point  which  appears  in  supersonic  separated 
and  reattaching  flow  analyses  and  has  been  shown  to 
be  analogus  to  the  critical  point  occurring  in  the 
throat  of  a subsonic-supersonic  nozzle.  Examples  of 
the  latter  type  are  the  "velocity  profile  critical 
points"  discussed  by  Shammroth  and  McDonald  [40-41]  and 
further  by  Gerhart  [34].  The  latter  type  of  critical 
point  is  identified  with  the  failure  of  a finite  set 
of  integral  equations  coupled  with  a particular  velocity 
profile  function  to  produce  a completely  independent 
set  of  equations  at  all  times.  Calculation  difficulties 
associated  with  the  latter  type  of  singularity  can  be 
overcome  by  using  an  over-constrained  set  of  equations 
[40,41]  or  by  switching  from  the  weak  to  strong  inter- 
action formulation  or  vice  versa  [34].  It  is  important 
that  these  latter  types  of  singularities  not  be  assigned 
physical  significance  and  it  is  questionable  if  smooth 
passage  through  such  singular  points  should  be  made  a 
criterion  for  arbitrary  adjustment  of  flow  parameters 
(see  Shammroth  and  McDonald's  [40,41]  discussion  of  the 
work  of  Green  [ 42]»  and  also  Tai  [43,44]). 

In  the  current  work,  both  types  on  singularity  were 
encountered.  As  pointed  out  by  Kuhn  and  Nielsen  [3], 
if  the  wall-wake  velocity  profile  (equation  3.5)  is 
assumed,  the  coefficients  of  dX/dx  in  all  of  the  differ- 
ential equations  (3.8-3.12)  vanish  when  X * 


* 0, 


which  becomes  a singularity  in  either  the  weak  or  strong 
interaction  formal ations.  This  behavior  may  be  identi- 
fied with  the  separation  point  singularity  and  has  some 
physical  implications.  Examination  of  the  wall-wake 
velocity  profile  function  shows  that  all  influence  of 
the  (zero)  skin  friction  disappears.  Downstream  of  the 
point  of  zero  skin  friction,  the  boundary  layer  may 
develop  as  a recirculating  flow  above  a solid  surface, 
in  which  case  the  skin  friction  attains  negative  values 
or  as  a wake  type  flow  in  which  the  skin  friction  re- 
mains zero;  therefore,  "critical"  point  behavior  is  not 
inconsistent. 

Whatever  its  cause,  if  a singularity  occurs  in  the 
calculations,  a method  must  be  devised  to  calculate 
through  it  or  jump  over  it.  Kuhn  and  Nielsen  [ 3 ] indi- 
cate that  the  singularity  can  be  removed  from  their 

formulation  which  uses  6*  rather  than  O as  the  strong 

★ 

interaction  parameter  . They  show  that  the  ratios  of 
coefficients  of  dX/dx  in  their  differential  equation 
set  are  finite  so  that,  by  dividing  the  momentum  and 
moment  of  momentum  by  the  differentiated  displacement 
thickness  definition,  two  independent  differential  equa- 
tions are  recovered.  These  are  solved  together  with  a 
prescribed  5*(x)  and  the  5*  definition.  Of  course,  using 


* Since  it  occurs  where  * 0 (fully  developed  separa 
tion),  this  singularity  must  always  be  dealt  with  in 
the  strong  interaction  formulation. 
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the  6*  definition  twice  does  not  introduce  any  new 
information,  it  simply  represents  an  artificial  way 
of  removing  a computational  difficulty.  Kuhn  himself 
states  that  the  method  is  essentially  equivalent  to 
prescribing  X(x)  as  the  interaction  parameter  as  in 
his  previous  work  [2].  Had  it  been  deemed  necessary, 
a similar  development  could  have  been  undertaken  in 
the  present  work,  however  extensive  calculations  re- 
vealed that  the  singularity  associated  with  the  vanish- 
ing of  \ is  very  weak.  A singularity  is  of  course 
associated  with  the  vanishing  of  the  determinant  of  the 

matrix  of  the  coefficients  of  the  derivative  dM  /dx, 

e 

dX/dx,  d6/dx,  du./dx  of  equations  3.8-3.12.  The  pre- 

P 

dictor-corrector  method  used  for  the  solution  of  the 
differential  equations  is  sensitive  to  the  approaching 
of  a singularity,  with  the  marching  step  size  being 
automatically  reduced  as  the  determinant  decreases.  Cal- 
culations have  shown  that  the  skin  friction  parameter 

X,  can  be  reduced  to  a value  of  t .0000005  without 

★ 

any  notable  effect  on  the  calculations.  The  deter- 
minant thus  behaves  as  shown  in  Figure  3.7,  with  no  zero 
crossing  by  the  determinant.  Since  this  singularity  is 
so  hard  to  detect,  it  is  extremely  unlikely  that  it 
would  ever  effect  the  calculations  significantly.  In 


* The  existence  of  the  singularity  was  nominally 

verified  by  setting  X»0.0,  in  the  computer  program 
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order  to  account  for  the  possibility  of  its  occurrence, 
it  was  decided  that  if  this  singularity  does  occur,  the 
values  of  the  derivatives  from  the  last  upstream  sta- 
tion are  used,  thus  the  variables  are  linearly  extrapo- 
lated through  the  singularity. 

Singularities  of  the  "non  significant"  type  also 
occur  in  the  current  method.  Typically  these  singular- 
ities are  not  associated  with  either  a row  or  column  of 
the  differential  equation  matrix  vanishing  but  simply 
the  appearance  of  a zero  determinant.  Singularities  of 
this  type  involve  a "zero  crossing"  determinant  as 
shown  in  Figure  3.8.  In  such  a case,  the  approaching 
singularity  is  detected  by  the  differential  equation 
solver,  with  the  result  that  calculations  usually  break- 
down before  the  singular  point  is  reached.  Numerical 
experimentation  has  revealed  the  following  about  these 
types  of  singularities: 

(1)  They  are  more  prone  to  occur  in  compressible 
flow  calculations,  becoming  stronger  and  more 
plentiful  as  the  Mach  number  is  Increased, 

(2)  The  singularities  only  occur  for  X negative 

or  very  small  positive  (i.e.  near  fully  devel- 
oped separati on ) , 

(3)  Although  these  singularities  occur  in  both 
the  weak  and  strong  interaction  formulations. 
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at  the  same  values  of  ©,  X,  6,  Ug)  iji 
both  formulations. 


The  latter  fact  provides  the  key  to  avoiding  these 
singularities;  all  that  is  necessary  is  to  switch  from 
weak  to  strong  interaction  calculations  (or  vice  versa) 
as  a singularity  is  approached.  This  behavior  and 
remedy  were  discussed  by  Gerhart  [34]  although  he  ap- 
parently had  not  uncovered  the  entire  picture. 

During  the  course  of  the  calculations,  the  deter- 
minant of  the  matrix  of  the  differential  equations  is 
monitored.  If  the  ratio  of  the  current  value  to  a 
reference  value  is  less  than  0.2,  it  is  assumed  that  a 
"zero  crossing"  or  singularity  is  eminent  and  the  calcu- 
lation is  switched  from  weak  to  strong  interaction  or 
vice  versa.  For  positive  values  of  X,  the  reference 
determinant  is  the  value  at  the  last  point  where  > .0015, 
for  negative  values  of  X,  it  is  the  initial  value  of 
the  determinant  at  the  last  mode  switch.  In  practice, 
the  choice  of  the  weak  or  strong  interaction  mode  is 
usually  dictated  by  the  criterion  of  equations  3.34  and 
3.35,  with  singularity  appearance  dictating  the  switch 
from  strong  to  weak  interaction  only  in  the  regions  of 
reverse  flow  (after  separation,  during  free  shear 
I layer  development). 
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At  this  point,  a further  discussion  of  the  method 
of  solving  the  differential  equations  of  the  viscous 
layer  is  in  order. 

These  equations  have  the  form 


dM^ 
‘xi  dT 


A + A ^ + 

''ai  dx  '^31  dx 


du 
i dlT 


6 


1 .4 


In  the  strong  interaction  formulation,  is  an  unknown 

and  there  are  4 differential  equations  to  solve;  in  the 

weak  interaction  formulation,  M is  known  and  since  O 

e 

does  not  appear  in  a derivative,  it  can  be  eliminated 
with  the  result  that  only  3 differential  equations  re- 
main. The  predictor-corrector  method  employed  for  the 
solution  of  the  differential  equations  requires  know- 
ledge of  the  derivatives  of  4 previous  points  to  march 
the  solution  forward;  therefore,  it  must  be  started  by 
a Runge-Kutta  method,  which  turns  out  to  be  a rather 
time  consuming  and  sensitive  process.  Now  if  in  switch- 
ing from  weak  to  strong  interaction,  we  add  an  extra 
differential  equation  to  the  set,  it  is  obvious  that 
the  calculations  must  be  restarted.  In  order  to  avoid 
this,  both  weak  and  strong  interaction  calculations  are 
arranged  to  solve  4 differential  equations;  when  making 
weak  interaction  calculations,  the  continuity  equation 
(i  »1  above)  is  subtracted  from  the  momentum  and  momentum 
equations  (i  *2,3)  and  discarded.  It  is  replaced  by 
the  equation 
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dM  . 

dx”  “ Tx  given  function] 

Thus  A =1,A  *A  ~ ^ =0,  B=  known 

1>1  lf2  1>3  li4  1 

function  of  x.  At  each  station,  after  a solution  has 
been  obtained,  the  criteria  for  switching  modes  are 
all  checked,  a decision  as  to  whether  to  advance  the 
calculations  via  the  weak  interaction  or  strong  inter- 
action formulation  is  made,  and  the  differential  equa- 
tions are  set  up  accordingly. 

3.7  Verification  of  the  Boundary  Layer  Method 

Since  exact  solutions  are  lacking,  the  ultimate 
test  of  any  (separating)  turbulent  flow  calculation  pro- 
cedure is  confrontation  with  experimental  data.  Since, 
for  attached  flows,  in  the  weak  interaction  mode,  the 
current  method  is  quite  similar  to  several  of  those 
presented  at  the  Stanford  Conference  [19,23,25],  it 
might  be  expected  to  perform  in  roughly  the  same  manner; 
preliminary  calculations  verified  this. 

In  order  to  verify  the  ability  of  the  method  to 
make  calculations  of  separating  boundary  layers,  it  was 
necessary  to  calculate  such  flows.  The  flows  to  be  cal- 
culated should  meet  the  following  criteria 

(A)  Both  incompressible  and  subsonic  (or  transonic) 
flows  should  be  considered. 


(B)  Extensive  details  of  the  measured  flow  are 
needed,  especially  including  both  pressure 
(Mg)  and  edge  angle  (0)  distributors,  as  well 
as  and  6 information, 

(C)  Both  two-dimensional  and  axisymmetric  geo- 
metries are  desirable, 

(D)  Both  separation  and  reattachment  of  the  experi- 
mental flow  are  desirable, 

(E)  Separation  must  be  caused  by  adverse  pressure 
gradients,  not  by  sudden  geometry  changes  such 
as  back  steps. 

Few  experimental  flows  following  satisfying  all  of 
these  criteria  were  found;  however,  two  excellent  test 
cases  were  found. 

Simpson  et  al.  [21]  have  made  exhausti ve  measure- 
ments in  an  incompressible,  two  dimensional  separating 
boundary  layer.  Then  experimental  results  include  both 
velocity  and  edge  angle  distributions,  boundary  layer 
thicknesses,  form  factors,  skin  friction  coefficients, 
and  details  of  the  turbulence  structure.  In  fact  the 
on^y  thing  which  keeps  this  flow  from  being  a complete 
incompressible  test  case  is  the  fact  that  the  flow  does 
not  reattach. 

Alber  et  al . [13]  have  presented  measurements  made 
in  a two-dimensional  transonic  separating  boundary  layer 
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flow.  They  present  information  on  pressure  distri- 
bution, boundary  layer  thicknesses,  form  factors  and 
skin  friction  coefficient.  Although  edge  angles  were 
not  presented  in  the  paper,  they  could  be  estimated 
from  given  6,  6*,  data  via: 

tan  © » ^ + (M*-l)(6-6  ) ;^  ^ (3.36) 

' 0 e 

This  flow  also  included  reattachment. 

The  two  flows  mentioned  satisfy  all  of  the  desirable 
criteria  except  for  axisymmetric  geometries.  Putnam 
and  Abeyounis  [45]  surveyed  the  flow  field  in  the  vicin- 
ity of  a boattailed  axisymmetric  afterbody  and  presented 
Mach  numbers  and  flow  angles;  however,  no  boundary  layer 
details  were  obtained  and  the  published  figures  giving 
M and  "©"  were  too  small  to  be  useful.  Attempts  to 
obtain  larger  scale  figures  were  unsuccessful. 

Using  the  method  developed  in  Sections  3.1  -3.5,  the  Simpson 
and  Alber  flows  were  calculated.  The  results  are  shown  in 
Figures  3.9-3.16.  The  calculations  were  started  by  match- 
ing Mg,  5,  X at  the  furthermost  upstream  point  shown. 

In  both  cases,  the  entire  experimentally  determined  Mg(x) 
and  0(x)  distributions  were  made  available  so  that  the 
parameter  required  by  the  choice  of  interaction  mode  was 
always  available.  The  normal  stress  terms  as  formulated 


* I.  £.  Alber  has  been  kind  enough  to  send  the  authors 
extensive  details  of  the  experimental  data. 


FIGURE  3.11  MEASURED  AND  PREDICTED  BOUNDARY  LAYER  THICKNESS 


DISTRIBUTION  FOR  SIMPSON  ET  AL.  FLOW 
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FIGURE  3.16  MEASURED  AND  PREDICTED  SHAPE  FACTOR 

DISTRIBUTION  FOR  THE  ALBER  ET  AL . FLOW 


FIGURE  3.13 


MEASURED  AND  PREDICTED  MACH  NUMBER  AND 
EDGE  ANGLE  DISTRIBUTIONS  FOR  THE  ALBER 
ET  AL.  FLOW 


FIGURE  3.14 

MEASURED  AND  PREDICTED  SKIN  FRICTION 

COEFFICIENT  DISTRIBUTION  FOR  THE 

ALBER  ET  AL.  FLOW 
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in  equation  3.33  ff  were  included  in  both  cases  and 
the  algebraic  turbulence  model  of  equations  3.24-3.26 
was  employed.  * Although  the  calculations  do  not  show 
perfect  agreement  the  following  are  believed  to  be  sig- 
nificant accomplishments. 

(1)  The  agreement  is  reasonable, 

(2)  Both  compressible  and  incompressible  flows 
are  predicted  by  a single  method, 

(3)  Separation,  reverse  flow,  negative  shear, 
and  reattachment  are  all  evident  and,  given 
the  difficulty  of  even  measuring  negative 
shear,  are  considered  rather  accurate, 

(4)  The  pressure  (Mg)  distribution  is  calculated 
reasonably  well  by  the  strong  interaction 
method . 

Based  on  the  results  of  these  trial  calculations  and 
in  comparison  with  the  experience  of  other  investigators 
in  calculating  these  flows  [27,33],  it  was  concluded 
that  the  integral  method  developed  herein  is  sufficiently 
accurate  to  be  used  as  part  of  an  a priori  separated 
flow  prediction  procedure. 


* It  should  be  pointed  out  that  slightly  better  calcu- 
lations could  be  obtained  by  adjusting  the  constants 
in  equation  3.24  for  each  flow.  These  "constants" 
might  be  a function  of  Mach  number  but  information  is 
insufficient  to  pursue  this  point  further. 


4.  Inviscid  Flow  Calculation  Procedure 

The  boundary  layer  method  described  in  the  previous 
chapter  represents  only  half  (or  less!)  of  the  entire  cal- 
culation procedure  for  a priori  prediction  of  separating 
boundary  layer  flows.  In  this  section,  results  of  efforts 
aimed  at  finding  solutions  to  the  inviscid  "half"  of  the 
flow  will  be  described. 

Methods  for  calculating  the  inviscid  (compressible  or 
incompressible)  flow  over  arbitrary  prescribed  plane  two- 
dimensional  and/or  axisymmetric  bodies  are  quite  plentiful 
and  need  not  be  described  here.  These  methods  basically 
solve  the  partial  differential  equations  2.4  and  2.5,  sub- 
ject to  the  boundary  conditions  of  2.7,  2.8  and  2.10.  The 
boundary  conditions  of  2.10  is  applied  at  the  surface  of 
the  body  in  question,  it  being  assumed  that  the  boundary  layer 
is  nonexistent,  in  which  case  © is  zero;  that  is  the  slope  of 
the  velocity  vector  is  required  to  be  equal  to  the  body  slope. 
This  constitutes  the  so  called  "direct"  or  "analysis"  prob- 
lem of  inviscid  flow  theory. 

A few  authors  have  presented  methods  for  calculating 
the  body  which  will  yield  an  arbitrarily  prescribed  pressure 
distribution  over  its  surface  [46-48].  Such  methods  formally 
solve  equations  2.4  and  2.5  subject  to  the  boundary  conditions 
of  2.7,  2.8  and  2.9  (again,  typically  6 = 0 so  the  boundary 

condition  is  applied  at  the  body  surface).  In  this  method, 

i 
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it  is  actually  the  slope  of  the  body  surface  which  is  computed, 
with  the  geometry  being  determined  by  an  integration  of  the 
surface  slope.  This  constitutes  the  "inverse"  or  "design" 
problem  of  inviscid  flow  theory. 

An  inviscid  flow  method  satisfactory  for  coupling  with 
the  separating  boundary  layer  method  of  Chapter  3 must 
combine  elements  of  both  the  direct  and  inverse  problems. 

We  have  seen  that,  at  most  stations  along  the  flow  surface, 
only  one  formulation,  weak  or  strong,  of  the  boundary  layer 
equations  will  apply.  Now  the  weak  interaction  formulation 
requires  that  be  specified  from  outside  information  (e.g. 
the  inviscid  flow  procedure)  and  computes  an  updated  velocity 
slope©,  for  handing  to  the  inviscid  flow  procedure;  on  the 
other  hand  the  strong  interaction  formulation  requires  that 
© be  specified  and  calculates  an  which  is  handed  to  the 
inviscid  flow  procedure.  It  is  therefore  obvious  that  the 
inviscid  flow  procedure  must  be  capable  of  accepting  mixed 
boundary  conditions,  with  direct  type  (9  specified)  condi- 
tions at  some  points  and  inverse  type  (fl^  specified)  condi- 
tions at  other  points  and  that  which  points  are  of  which  type 
are  predetermined  by  the  boundary  layer  calculations.  Two 
other  complications  arise.  First,  because  the  boundary  layer 
may  change  from  iteration  to  iteration,  the  type  of  boundary 
condition  to  be  applied  at  a particular  point  may  change  from 
iteration  to  iteration.  Second,  the  boundary  conditions  are 
not  applied  on  a "solid  surface"  which  is  both  impermeable 
and  fixed  from  iteration  to  iteration  but  instead  on  a surface* 

* Hereafter  called  the  "6  surface" 
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representing  the  edge  of  the  boundary  layer  (the  surface  j 

! 

3 

y » 6(x),  r ■ R + 5(x)  cos  a(x))  which  Is  neither  Impermeable  i 

(not  a significant  problem)  nor  fixed  from  Iteration  to  ^ 

iteration  (a  significant  problem). 

It  should  not  be  surprising  that  no  currently  available 
Inviscid  flow  procedure  Incorporating  all  of  these  require- 
ments was  found;  therefore.  It  was  necessary  to  attempt  to  | 

develop  a satisfactory  procedure.  In  the  research  efforts  \ 

described  herein,  3 different  methods  were  Investigated.  | 

These  were 

(1)  The  method  of  Integral  relations 

(2)  The  surface  singularity  method 

(3)  The  finite  element  method 

The  formulations  for  the  methods  of  Integral  relations 
and  finite  elements  were  done  for  either  plane  two-dimensional 
or  axisymmetric  bodies.  The  formulation  of  the  surface  source 
method  was  for  axisymmetric  geometries  only. 

Each  of  these  methods  will  be  the  subject  of  a separate 
section. 

4.1  The  Method  of  Integral  Relations 

The  method  of  Integral  relations  was  the  method 
Initially  Investigated  for  the  solution  of  the  external  compres- 
sible flow  problem.  As  the  method  Is  conceptually  and 
numerically  similar  to  integral  boundary  layer  methods, 
it  appeared  to  have  certain  advantages  over  more  conven- 
tional Inviscid  flow  methods  for  the  interaction  problem. 

■ 
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Briefly,  the  method  involves  Integration  of  the 
system  of  flow  equations  in  the  transverse  direction. 

To  perform  this  integration,  the  transverse  variation 
of  the  integrands  must  be  known.  A general  approach 
is  to  approximate  this  variation  with  polynomials. 

The  result  is  that  the  original  elliptic  partial  dif- 
ferential equations  of  the  flow  are  reduced  to  a para- 
bolic two-point  boundary  value  problem  that  may  be 
solved  numerically  using  a standard  Runge-Kutta  or 
predictor-corrector  subroutine,  also  needed  for  the 
viscous  flow  problem. 

It  was  hoped  that  reduction  of  the  problem  to 
ordinary  differential  equations  would  result  in  solu- 
tion scheme  that  was  consistent  with  the  integral 
boundary  layer  method,  able  to  handle  arbitrary  flow 
geometries,  and  still  be  simpler  and  faster  than  finite 
difference  methods. 

Experience  with  the  method  of  integral  relations 
eventually  demonstrated  that  any  advantages  gained  by 
reduction  of  the  problem  to  ordinary  differential  equa- 
tions were  outweighed  by  the  algebraic  complexity  of  the 
resulting  system.  The  equations  were  difficult  to 
derive  and  program,  and  ran  slowly.  Furthermore,  the 
iteration  schemes  required  to  handle  the  two-point 
boundary  valve  problems  were  physically  unrealistic. 

For  this  and  other  reasons  detailed  in  the  following. 
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research  on  the  method  of  integral  relations  was  dis- 
continued in  favor  of  the  finite  element  method. 

4.1.1  Literature  Review 

A.  A.  Dorodnicyn  first  introduced  the  method  of 
integral  relations  in  1959  for  mixed  elliptic- 
hyperbolic  aerodynamic  problems  [49].  He  and  his 
colleagues  solved  subsonic  flow  over  ellipses  and  ellip- 
soids, transonic  flow  over  an  ellipse,  and  supersonic 
flow  over  a cylinder,  although  his  paper  presents  few 
details.  Holt  of  the  University  of  California  at 
Berkeley  used  the  method  to  solve  the  transonic  flow 
over  a cylinder  [50].  Melnik  and  Ives  of  Grumman  Aero- 
space solved  compressible  flows  over  a cylinder,  an 
ellipse,  and  simple  non-lifting  airfoil  sections  using 
the  method  of  integral  relations  [51]. 

T.  C.  Tai  of  the  Naval  Ship  Research  and  Develop- 
ment Center  has  used  the  method  for  more  practical  prob- 
lems. He  reports  considerable  success  using  the  method 
to  solve  the  supercritical  flow  over  symmetric  airfoils 
[5]]  and  lifting  airfoils  [53],  and  matched  to  laminar 
and  turbulent  boundary  layer  computations  [43-44].  Tai's 
work  is  impressive,  but  a typical  airfoil  problem  takes 
considerable  computer  time  and  requires  interactive 
graphics  capabilities  [54]. 
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4.1.2  Formul at  ion 

The  basic  equations  for  an  inviscid  compres- 
sible fluid  in  Cartesian  or  cylindrical  coordinates 
can  all  be  written  in  divergence  form  as  follows: 

Continuity 

ll  '“''z*  I?  (“''r*  * J V " 

r-momentum 

It  1'^'’  * »''!>  * If  '^’'zV*  - io  ^ ■ 0 

r-momentum 

It  + If  (KP  + pvp  + j — - 0 

K = 1/yM^  ® constant 
00 

All  velocities  are  non-dimensional i zed  by  lengths 
by  some  typical  length  (usually  the  body  length),  and 
pressure  and  density  by  their  free-stream  values. 

Further,  density  may  be  related  to  the  velocity 
field  by  the  energy  equation,  and  pressure  may  be  related 
to  density  by  the  isentropic  relation.  For  incompressible 
flow,  0=1,  and  pressure  is  related  to  velocity  by  the 
Bernoulli  equation. 

Finally,  the  boundary  conditions  are  as  follows: 

At  a solid  surface,  the  normal  velocity  is  zero.  At 
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"infinity",  the  flow  is  undisturbed,  i.e.  * p = 

P = 1,  * 0. 

4.1.3  Development  of  the  Method  of  Integral  Relations 
Equations 


To  apply  the  method  of  integral  relations,  the 
flow  equations  must  be  written  in  divergence  form: 

Tz  A(^,r,V2.  B(z,r,v^. ..  ) = R(z,r,v^...) 


Like  an  integral  boundary  layer  method,  the  method  of 
integral  relations  relies  on  integration  of  the  flow 
equations  in  the  r direction,  reducing  the  partial  dif- 
ferential equations  to  ordinary  ones  with  independent 
variable  z.  To  do  so,  the  r-variation  of  the  inte- 
grands must  be  known.  If  the  flow  field  is  divided 
into  strips  bounded  by  typical  streamlines,  the  inte- 
grands may  be  approximated  by  polynomials  of  the  form: 

n . N . 

A = I a.(r-r  , R = S b.(r-r  )’ 
i*0  ^ ° i=0 


where  N is  the  number  of  strips  and  a^  and  b are  con- 
stants evaluated  on  strip  boundaries. 

After  tedious  integration  and  rearrangement,  the 
equations  can  be  reduced  to  the  forms: 

3T-  ■ dT"  • dz"  • ' ' > 


81 


V 3 W ■ ■ !■ 

ro  zo  dz 


^ . f (V  V iii  1 
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a a (u  V clAn-1  x 
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where  the  subscript  {)^  refers  to  the  stagnation  stream- 
line, and  the  subscript  ()^  is  the  number  of  a strip 
boundary  as  shown  in  Figure  4.1.  The  full  equations  for 
two-dimensional  flow  may  be  found  in  [53]]. 

It  is  unreasonable  to  expect  to  approximate  all 
flow  quantities  from  r^  to  r^  with  a simple  polynomial. 
Instead,  the  flow  domain  is  treated  as  a series  of  ef- 
fective strips,  using  a second-order  polynomial  approx- 
imation across  each  pair  of  strips,  as  shown  in  Figure 
4.1.  Given  a pair  of  strips  bounded  by  r^ , r^_^  , and  r^, 
the  flow  equations  can  be  integrated  along  r^_i  if  the 
flow  properties  are  known  at  r^.  Thus,  the  solution 
can  be  found  along  r^  given  properties  on  r2-  Similarly, 
the  solution  can  be  found  along  rg  given  properties  on 
r^,  and  so  on  until  the  free  stream  boundary  is  reached 
where  the  flow  properties  are  known.  The  integration  of 
all  sets  of  equations  is  carried  out  simultaneously  at 
successive  z-stations.  Strips  may  be  added  at  any  point 
in  the  flow  to  increase  accuracy. 


0 
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FIGURE  4.1  A MULTI-STRIP  SCHEME  IS  USED  TO  APPLY  THE  FREESTREAM 
BOUNDARY  CONDITIONS  AND  TO  IMPROVE  ACCURACY.  THE 
METHOD  OF  INTEGRAL  RELATIONS  REDUCES  THE  PARTIAL 
DIFFERENTIAL  EQUATIONS  OF  THE  FLOW  TO  ORDINARY  ONES 
ALONG  STREAMLINES.  FREESTREAM  CONDITIONS  ACT  AS 
BOUNDARY  CONDITIONS  FOR  THE  EQUATIONS  ALONG  . 

THE  SOLUTION  ALONG  THEN  BECOMES  BOUNDARY  CONDITIONS 
FOR  THE  EQUATIONS  ALONG  r^,  ETC. 
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4.1.4  Iteration  Procedures  for  Subsonic  Flow  Over 
Symmetric  Bodies 


Integration  reduces  the  elliptic  partial  dif- 
ferential equations  of  the  flow  to  ordinary  differen- 
tial equations.  Two  iteration  schemes  used  to  solve 
these  boundary  value  problems  preserve  the  elliptic 
nature  of  the  flow. 

Flow  integration  starts  well  upstream  where  free- 
stream  conditions  are  imposed.  A mathematical  dis- 
turbance, based  on  predicted  stagnation  conditions, 
must  be  applied  along  the  stagnation  streamline  for  the 
flow  to  vary  at  all.  Details  may  be  found  in  [53]. 

When  the  flow  begins  to  vary  stably,  the  disturbance 
is  removed,  and  integration  proceeds  to  the  body. 

There,  the  predicted  conditions  may  be  evaluated,  and 
the  new  values  are  used  to  start  another  iteration. 

This  procedure  repeats  until  convergence  is  achieved. 

Integration  proceeds  from  stagnation,  along  the 
body,  and  well  downstream,  where  the  pressure  is  expected 
to  return  to  its  free-stream  value.  It  turns  out  that 
this  condition  is  extremely  dependent  on  the  location 
of  the  free-stream  boundary  r^,  apparently  since  the 
polynomial  approximations  involve  r^.  Hence,  iterating 
on  r^  to  force  the  downstream  pressure  to  converge  pro- 
vides feedback  from  downstream  to  upstream. 
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4.1.5  Results  and  Conclusion 

The  upstream  iteration  for  a symmetric  Joukowsky 
airfoil  at  * 0.5  was  completed  on  a simple  two  strip 
grid.  In  general,  the  flow  behavior  was  correct,  but 
cannot  be  compared  to  other  solutions  since  the  down> 
stream  iteration  never  converged  to  an  "exact"  value 
of  r^.  Similar  test  runs  were  made  using  multiple  strip 
grids.  From  these  tests,  the  computer  time  required  for 
the  solution  of  the  flow  over  a typical  b'^dy  was  pre- 
dicted to  be  seven  minutes  on  an  IBM  370/158. 

Several  objections  to  the  method  of  integral  rela- 
tions led  to  termination  of  the  research  before  any 
complete  solutions  were  obtained. 

An  initial  objection  was  that  the  integrated  equa- 
tions are  algebraically  cumbersome.  Second  order  approxi- 
mations are  very  difficult  to  derive  and  program,  third 
order  would  be  prohibitive. 

A second  objection  is  that  although  the  equations 
are  in  primitive  variable  form,  that  is,  the  unknown 
quantities  are  the  velocities,  density,  and  pressure 
rather  than  a potential  function,  the  quantities  solved 
for  are  combinations  of  these  variables  and  must  be  de- 
coupled algebraically  at  each  step. 

The  major  objection  is  that  the  iterative  schemes 
required  to  preserve  the  elliptic  nature  of  the  flow 
are  artificial.  Upstream  it  is  annoying  that  a disturbance 
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Is  required  to  cause  the  flow  to  vary  at  all.  Near 
the  body,  the  method  of  Integral  relations  equations 
become  unstable,  and  the  solution  must  be  extrapolated 
to  the  stagnation  point.  It  is  comforting  to  note  that 
the  upstream  disturbance  is  based  on  physically  realistic 
stagnation  properties,  and  that  the  upstream  iteration 
converges  quickly.  Regarding  the  downstream  iteration, 
it  is  unrealistic  that  the  location  of  the  free  stream 
boundary  should  have  a pronounced  effect  on  convergence 
of  the  downstream  pressure,  yet  this  is  a most  sensitive 
iteration. 

Finally,  the  large  computing  times  predicted  for  a 
simple  airfoil  demonstrated  that  the  method  of  integral 
relations  is  not  competitive  with  finite  difference  or 
finite  element  methods. 

4.2  Surface  Singularity  Method 

When  the  research  described  herein  was  begun,  it 
was  decided  that  of  all  inviscid  flow  methods  available, 
the  surface  singularif'  'ethod  held  the  most  promise  of 
being  easily  adopted  to  the  special  needs  of  separated 
flow  computation.  This  method  has  undergone  extensive 
development  by  researchers  at  the  McDonnel 1 -Douglas 
Corporation  [55-58].  The  work  by  Hess  and  Smith  [55] 
may  be  considered  classic.  Computer  programs  for  solving 
the  direct  problem  are  widely  available  and  modifica- 
tions to  solve  the  inverse  problem  have  been  documented 
[48,52,58]. 
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Basically,  the  surface  singularity  method  constructs 
a solution  to  the  velocity  potential  differential  equa- 
tion, which  replaces  2.4  and  2.5  by  the  single  equation 
(the  adiabatic-isentropic  energy  equation  is  also 
i ncl uded ) 


The  velocity  isv  = V cos  0,  v^  = t^*V  sin  9 

z oz  r or 

This  equation  is  equivalent  to  2.4  and  2.5 
for  incompressible  flow  (M^=0)  and  a small  perturba- 
tion approximation  for  > 0.  The  equation  is  apparently 
accurate  up  to  M = 0.8. 

It  can  be  shown  that  4.1  can  be  satisfied  by  various 
singularities  such  as  sources,  sinks,  doublets,  and  for 
plane  flow,  vorticies.  Since  the  equation  is  linear, 
combinations  or  distributions  of  singularities  also 
satisfy  the  equation.  In  the  surface  singularity  method, 
singularities  are  distributed  over  surfaces  (usually 
corresponding  to  body  surfaces)  in  the  flow  field.  The 
singularities  generate  a flow  field  which  satisfies 
(4.1)  identically;  the  strengths  of  the  singularities 
(the  intensity  of  the  distribution)  are  determined  by 
requiring  that  the  flow  field  generated  satisfy  the 
boundary  conditions  (usually  tangency  at  a solid  surface). 
The  result  is  a Fredholm  integral  equation  of  the  second 
kind  for  the  unknown  singularity  density  function.  As 


7*“ 


87 


developed  by  Hess  and  Smith  [55],  the  corresponding 
numerical  procedure  involves  distributing  source 
(sink)  singularities  in  a stepwise  fashion  over  a 
piecewise  linear  approximation  (inscribed)  to  the 
actual  surface.  This  results  in  a set  of  linear  alge- 
braic equations  with  the  (piecewise)  source  densities 
ad  unknowns  as  an  approximation  to  the  integral  equa- 
tion. Subsequent  improvements  have  introduced  curved 
surface  elements,  doublet  or  combination  singularities, 
and  polynomial  singularity  distributions  over  each 
element  [56-58]. 

To  apply  the  method  to  the  current  problem,  we 
recognize  that  our  goal  is  to  calculate  the  flow  on 
and  external  to  a predetermined  surface  in  the  flow 
field,  the  surface  describing  the  edge  of  the  boundary 
layer.  This  surface  is  given  by 

r(z)  » R(z)  + 6(z)  cos  a(z) 

External  to  this  surface,  it  is  assumed  that  the  flow  is 
irrotational  and  inviscid,  thus  the  flow  satisfies  4.1. 


Introducing  the  transformations 


Z » / l -M*  z„  » 8z. 

oo  0 0 


4i  * 1 /S  41, 


results  in 


_£  + i_ 

3rJ  ^ 


the  "incompressible"  form  of  the  equation.  The 
"incompressible"  velocities  are  related  to  the  compres- 
sible velocities  by 


V +V 
zo  0 


1 - 6^  + ^ cos  9 


6 ^ sin  e 
0 00 


The  "6  surface"  is  first  modified  by  replacing  x by  , 
transforming  to  an  "incompressible"  geometry.  Next 
the  "5  surface"  is  divided  into  n segments  and  (source) 
singularities  are  distributed  in  stepwise  fashion  on 
the  "5  surface".  Once  these  singularities  have  been 
distributed,  we  can  form  the  influence  coefficient 
matrix  following  identically  the  steps  and  calculations 
of  Hess  and  Smith  [55].  In  the  current  method,  computer 
subroutines  adapted  directly  from  a computer  program 
written  at  the  Douglass  Aircraft  Company  ("EODA")  were 
used.  The  influence  coefficient  matrices  are  written 
and  Y^j  and  are  the  z and  r components  of  velocity 
induced  at  the  midpoint  of  the  "i"th  surface  element  by 
a unit  source  density  on  the  "j"th  element.  If  a.  is 

J 

the  actual  source  density  on  the  "j"th  element,  then 
the  total  induced  velocities  on  the  "i"th  element  are 
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Introducing  the  relationship  between  the  "incompres- 
sible" and  actual  (compressible)  velocities: 

Z X.  .a.  = cos  6)  .-1  ] (4.2) 

1=1  ' j j 'oo  ' 


^ V 

z Y..a.  = 6 sin  6).  (4.3) 

i=l  '''  J ''oo  ^ 


Now  ^)  is  the  ratio  of  the  actual  velocity  at  point  i 

’'oo  i 

on  the  "5  surface"  to  the  actual  free  stream  velocity 

★ 

and  is  hence  related  to  M and  M , 9,.  is  the  angle  of 

the  velocity  vector  with  respect  to  the  axis  at  point  i 
on  the  "6  surface".  With  reference  to  Figure  4.2,  note 
that  8 is  not  necessarily  equal  to  ei ther  the  body  angle 
a o_r  the  boundary  layer  flow  angle  Q but  is  related  to 
them  by: 

9 = a + 0 (4.4) 

Note  that  is  there  were  no  boundary  layer,  the  "5 
surface"  would  correspond  to  the  body  surface,  0 would 
be  zero,  and  9 would  equal  the  body  angle  a. 

Now  consider  equations  4.2  and  4.3.  From  a previous 
boundary  layer  calculation,  <5(z))  is  known,  thus 

the  geometry  of  the  "6  surface"  can  be  determined  and 
the  influence  coefficient  matrices  X and  Y can  be  calcu- 
lated. At  each  point  on  the  "6  surface",  either 
(which  can  be  used  together  with  to  find  V/V^)  ojr  0 


♦ 


v/v. 


00 


1*  ^ "i 


1/2 
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(which  determines  0 via  4.4)  has  been  calculated  by 

the  boundary  layer  method.  Equations  4.2  and  4.3  thus 

provide  2n  equations  for  no's  (irrelevant  to  the  overall 

method)  and  n of  el  ther  V/V^  or_  9,  which  provide  new 

estimates  of  or  © for  use  in  a further  cycle  of 
e 

boundary  layer  calculations.  At  any  point,  it  is 
(theoretically)  possible  to  calculate  either  or  ©, 
given  the  other.  Note  that  the  variable  which  must  be 
calculated  at  a particular  point  might  change  with  each 
(overal 1 ) i teration. 

All  that  remains  is  to  construct  a method  for 

solving  4.2  and  4.3  which  allows  for  either  V/V^  or  6 

to  be  determined.  Now  if  all  of  the  9's  are  known, 

with  all  of  the  V/V  's  to  be  determined,  the  equations 

00 

are  1 i near ; in  addition,  if  (S  cos  9)  x (equation  4.3) 
is  subtracted  from  (sin  9)  x (equation  4.2),  the  result- 
ing equation  is  easily  solved  for  the  (n)a's;  the  velocity 
is  then  easily  calculated.  This  is  in  effect  identical 
to  the  applications  of  the  method  to  the  calculation  of 
the  flow  over  specified  arbitrary  geometries  as  developed 
by  Hess  and  Smith. 

If  any  of  the  0's  are  unknown  (a  case  which  arises 
at  any  point  at  which  we  were  obliged  to  use  the  strong 
interaction  boundary  layer  formulation),  2 non-linear 
equations  are  introduced  at  that  point  and  standard 
matrix  methods  can  no  longer  be  used.  In  a typical 


separated  flow  calculation,  6 will  have  to  be  calcu- 
lated from  the  Inviscid  flow  equation  at  several  points 
so  a method  of  solving  the  mixed  set  of  2n  linear  and 
non-linear  algebraic  equations  arising  from  4.2  and 
4.3  with  V/V  known  and  9 unknown  was  developed. 

Now  4.2  and  4.3  can  be  written: 


/1-S|  . 1 , 


(Xj)  . 0 


(4-5) 


where  V.  » V/V  ). 

S^.  = sin  9^ 

Xj  . Oj  j.1*n 

X.  * V.  or.  S.  j.*n+l-^2n 

j I I 

It  was  found  necessary  to  use  sin  9 as  the  unknown 
rather  than  cos  9 The  two  are  related  via  the  square 
root  identity  e.g.  s1n  9 =-/l  - cos^  e or 


cos  9 


- s i n2  9 Cos  9 can  be  assumed  to  always 


have  the  positive  sign  since  the  external  flow  Is  not 
expected  to  flow  "upstream";  however,  the  sign  of  sin  9 
cannot  be  determined  since  the  flow  could  have  an  "up" 
or  "down"  component.  Viewed  another  way,  although  we 
can  always  preassign  the  proper  sign  to  cos  9,  we  can- 
not to  sin  9 and  must  rely  on  the  calculations  to  set 
the  sign  of  sin  9. 


w 
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The  equations  (4.5)  are  solved  via  Newton 
Iteration: 


- x"!*  + Sx^ 
J 3 J 

where 


3F, 

3x 


-I 

j x*x 


and 


3F. 
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(4.6) 


(4.7) 


(4.8) 


If  Is  to  be  determined,  the  upper  term  appears  on 
the  right  half  of  the  matrix.  If  Is  to  be  determii 
the  lower  term  appears.  Solution  of  these  equations 


proceeds  as  follows:  First  5a.,  1s  solved  from: 

3 

(a  set  of  nxn  linear  algebraic  equations) 
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where  is  the  (nxn)  diagonal  matrix  in  the  upper 

right  hand  corner  of  4.8  and  is  the  (nxn)  diagonal 

matrix  in  the  lower  right  hand  corner  of  4.8.  Then 
6V^  or  5S^  (whichever  is  required)  is  calculated  from 

«v,  . (-f(')(x^)  - L.  EX,j4Cj)/a(',> 
or 

6S,  ■ (-F{^>(xp  - I V,j60j)/A(f 

It  was  found  necessary  to  calculate  5V^.  and  6S. 
using  the  separate  formulations  above.  This  is  felt 
to  be  due  to  the  "dominance"  of  the  "z " direction 
over  the  "r"  direction  in  the  velocity. 

The  development  outlined  here  is  equally  applic- 
able to  either  plane  2-dimensional  (non-lifting)  or 
axisymmetric  geometries.  The  only  difference  would 
appear  in  the  formulation  of  the  influence  coefficient 
matrices  X..  and  Y...  In  this  work,  since  the  senior 

I J I J 

author  was  considerably  more  familiar  with  the  axi- 
symmetric method,  only  axisymmetric  formulations  were 
programmed . 

The  "di rect- i nverse"  inviscid  flow  method  out- 
lined above  was  tested  in  the  following  manner.  An 
axisymmetric  body  geometry,  typical  of  those  used  in 
the  experimental  investigation  of  nozzle  afterbody  drag 
[45],  was  specified  (see  Figure  4.3).  The  inviscid 
flow  over  the  body  surface  was  calculated  by  specifying 
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the  angle  everywhere  and  calculating  the  velocities. 
(Equivalent  to  a standard  "direct"  problem.)  A number 
of  points  In  the  boattall  region  were  then  selected  as 
"Inverse  points"  and  a calculation  was  done  In  which 
the  velocities  at  these  points  were  assigned  the  values 
calculated  in  the  "direct"  calculation  with  the  angles 
to  be  determined,  while  at  the  remaining  body  points 
the  angle  was  again  specified.  The  initial  guesses 
were 

c,  . 0 

9^  = 0 at  points  where  9 is  to  be  determined 
V/V^)^  * 1 at  points  where  V is  to  be  determined 
i.e.  uniform  parallel  flow. 

The  computation  converged  to  the  correct  velocity 
and  geometry  in  7 Iterations.*  The  iteration  history 
is  shown  in  Figure  4.4.  The  (possibly  limited)  ability 
to  make  mixed  duct-inverse  inviscid  flow  calculations 
using  the  surface  source  formulation  was  believed  to  be 
verified  by  this  calculation. 

4.3  The  Finite  Element  Method 

The  finite  element  method  is  a numerical  technique 

that  originated  in  structural  analysis,  but  is  proving 

to  be  a powerful  tool  in  all  continuum  problems.  In 

the  finite  element  method,  the  domain  of  interest  is 

divided  into  many  smaller  domains  or  finite  elements. 

* The  values  of  V/V  and  a were  obtained  in  3 iterations,  the  values 

00 

of  9 converged  most  slowly. 
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Corners  of  the  elements  are  known  as  nodes.  The  de- 
pendent variables  of  the  problem  are  approximated  by 
interpolation  functions  across  each  element.  The 
Galerkin  method,  a subclass  of  the  method  of  weighted 
residuals,  is  used  to  minimize  the  error  resulting 
from  use  of  the  interpolation  functions  in  the  govern- 
ing equations.  Assembly  of  the  Galerkin  equations 
from  each  element  results  in  a global  system  of  algebraic 
equations  for  the  nodal  values  of  the  dependent  vari- 
ables, which  may  be  solved  by  standard  matrix  methods. 

A finite  element  program  was  developed  to  solve 
inviscid  compressible  flows  over  arbitrary  two-dimensional 
or  axisymmetric  bodies.  Finite  difference  methods  are 


by  far  more  popular  for  this  type  of  problem,  but  the 
finite  element  method  has  certain  advantages.  The  most 
obv.uus  is  the  ability  of  the  method  to  fit  arbitrary 
geometries  by  judicious  placement  of  the  elements,  a 
necessity  in  the  boundary  layer  interaction  problem. 
Secondly,  work  by  Popinsky  and  Baker  [59]  indicates 
that  on  coarse  grids  finite  element  methods  are  more 
accurate  than  finite  difference  methods.  Third,  the 
inviscid  compressible  flow  problem  is  best  formulated 
in  terms  of  the  potential  function  for  finite  differ- 
ence solution.  As  will  be  discussed  later,  the  finite 


element  solution  is  more  conveniently  formulated  in 
terms  of  primitive  variables.  Two  unknowns,  the  velo- 
city components,  are  solved  for  directly  at  each  node. 


FIGURE  4.4  ITERATION  HISTORY  FOR  INVERSE  INVISCID  FLOW  CALCULATION  OF  BOATTAILED  BODY 
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Corners  of  the  elements  are  known  as  nodes.  The  de- 
pendent variables  of  the  problem  are  approximated  by 
interpolation  functions  across  each  element.  The 
Galerkin  method,  a subclass  of  the  method  of  weighted 
residuals,  is  used  to  minimize  the  error  resulting 
from  use  of  the  interpolation  functions  in  the  govern- 
ing equations.  Assembly  of  the  Galerkin  equations 
from  each  element  results  in  a global  system  of  algebraic 
equations  for  the  nodal  values  of  the  dependent  vari- 
ables, which  may  be  solved  by  standard  matrix  methods. 

A finite  element  program  was  developed  to  solve 
inviscid  compressible  flows  over  arbitrary  two-dimensional 
or  axisymmetric  bodies.  Finite  difference  methods  are 
by  far  more  popular  for  this  type  of  problem,  but  the 
finite  element  method  has  certain  advantages.  The  most 
obvious  is  the  ability  of  the  method  to  fit  arbitrary 
geometries  by  judicious  placement  of  the  elements,  a 
necessity  in  the  boundary  layer  interaction  problem. 
Secondly,  work  by  Popinsky  and  Baker  [59]  indicates 
that  on  coarse  grids  finite  element  methods  are  more 
accurate  than  finite  difference  methods.  Third,  the 
inviscid  compressible  flow  problem  is  best  formulated 
in  terms  of  the  potential  function  for  finite  differ- 
ence solution.  As  will  be  discussed  later,  the  finite 
element  solution  is  more  conveniently  formulated  in 
terms  of  primitive  variables.  Two  unknowns,  the  velo- 
city components,  are  solved  for  directly  at  each  node. 
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and  do  not  have  to  be  computed  from  derivitives  of 
a velocity  potential.  Finally,  the  finite  element 
method  tends  to  produce  "neat"  algorithms  that  are 
easy  to  program. 

The  computer  program  was  tested  extensively  against 

analytic  incompressible  solutions.  Excellent  results 

( 

were  obtained  for  the  pressure  coefficients  on  Rankine 
ovals  and  ovoids,  on  a sohere  and  a cylinder,  and  on  a 
14%  thick  Joukowsky  airfoil.  Compressible  results 
agreed  well  with  those  predicted  by  compressibility 
transformations.  The  pressure  coefficient  over  a NASA 
boattail  model  was  computed,  and  the  results  agreed  well 
with  a published  finite  difference  solution  [60]  over  a 
range  of  subsonic  Mach  numbers. 

The  inviscid  flow  program  was  also  coupled  iter- 
atively with  a Sasman-Cresc i integral  boundary  layer  pro- 
gram [61],  using  the  classical  method  of  augmenting  the 
body  by  the  displacement  thickness.  Results  for  the 
NASA  boattail  model  agreed  reasonably  well  with  pub- 
lished data  [60].  This  work  clearly  demonstrated  the 
ease  with  which  the  finite  element  method  can  be  made 
to  follow  a variable  geometry. 

Finally,  inverse  or  design  calculations  have  been  attempted 
in  which  the  inviscid  flow  program  was  to  be  used  iter- 
atively to  compute  the  axisymmetric  geometry  correspond- 
ing to  a prescribed  pressure  distribution.  Two  formu- 


lations  of  the  problem  were  tested,  the  sole  difference 
being  the  boundary  condition  applied  on  the  non- 
converged  body.  Only  one  method  gave  promising  re- 
sults, but  a converged  solution  was  not  obtained. 
Development  of  this  approach  is  continuing. 

4.3.1  Literature  Review 

Incompressible  ideal  flow  has  been  the  topic  of 
many  finite  element  papers.  Habashi  [10]  solved  lifting 
airfoil  problems  by  mapping  the  airfoils  to  near  circles, 
then  discreet! zing  the  resulting  finite  field  with  trian- 
gular elements  spanned  by  linear  interpolation  functions. 
As  a free-stream  boundary  condition,  Habashi  applied  the 
asymptotic  form  of  the  analytic  solution  for  flow  over 
a cylinder.  The  circulation,  and  hence  the  lift,  was 
solved  for  directly  as  a problem  unknown.  Habashi's 
program  is  efficient  and  accurate,  and  demonstrates  the 
utility  of  the  finite  element  method. 

T.  J.  Chung's  notes  from  the  short  course  "Finite 
Element  Methods  in  Fluid  Dynamics",  University  of  Alabama, 
Huntsville,  1976  [62],  include  complete  formulations  and 
computer  programs  for  two-dimensional  ideal  flow  using 
triangular  elements,  and  for  axisymmetric  flow  using 
quadrilateral  elements.  Other  papers  on  simple  ideal 
flows  include  those  by  Norie  and  de  Vries  [63],  Schmidt 
[64],  Shen  [65]  and  Street  [66]. 


An  excellent  paper  using  the  finite  element  method 
for  Inviscid  compressible  flow  Is  Hirsch's  computation 
of  turbomachine  through  flow  [67]. 

Chung  [62]  and  Heubner  [68]  have  both  presented 
finite  element  formulations  of  the  full  potential  equa- 
tion. Chung  and  Hooks  [69]  have  used  this  type  of 
formulation  to  obtain  some  Initial  shockless  results  for 
flow  over  a small  bump,  but  their  primary  objective  Is 
to  use  an  element  with  discontinuous  Interpolation 
functions  to  solve  flows  with  shock  waves.  Chung  and 
Chlou  have  formulated  unsteady  compressible  flow  In 
terms  of  the  equations  of  continuity,  momentum,  energy, 

; and  state,  using  primitive  variabl es  [20] . Using  this 

formulation,  they  have  solved  for  two  points  In  the 
unsteady  boundary  layer  behind  a moving  shock  wave,  with 
reasonable  results. 

The  most  common  use  of  the  finite  element  method 
In  compressible  flow  so  far  has  been  In  the  solution  of 
small  perturbation  forms  of  the  potential  equation. 

Carey  solved  the  Incompressible  flow  over  a cylinder 
then  used  this  result  to  obtain  a first  order  correction 
; [71,72].  Leonard  used  a similar  method  to  solve  the 

supersonic  flow  over  a Prandc l -Meyer  expansion  corner 
I [73].  Habashi  used  the  Prandtl -G1 auert  similarity 

form  of  the  potention  equation  to  solve  the  flow  over 
a cyl  1 nder  [74] . 

M 
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The  preceeding  references  Indicate  that  the  finite 
element  has  been  applied  at  least  Initially  to  many 
compressible  flow  problems.  However,  by  no  means  can 
It  be  said  that  the  method  has  been  fully  Investigated. 


~ woi'  <1  :_t2:^ 


r-’iT'h^  iti-'* 
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3P  , ^ ^ , gz  3P 

3x^  dp  3x^  3x^ 

along  with  the  equations  above,  to  reduce  the  z momentum 
equation  to  the  form 


!li  + !L:+  jL:.  V*  > I vM  ^ 

3z  3r  Jp  2 2 ''z  2 ''r*'  3z 


. vj] 


^^(v^  + v5-l)j^^2v,v 


z r 


3v. 


z r 3z 


} » 0 


(4.13) 

All  velocities  have  been  normalized  by  and  all  lengths 
by  some  arbitrary  length  Z. 

This  equation  is  to  be  solved  along  with  the  normal- 
ized irrotational ity  condition: 


3v^  ^ 

3r  3z“ 


(4.14) 


Equations  (4.13)  and  (4.14)  are  a coupled  set  of 

first  order  non-linear  partial  differential  equations 

for  the  normalized  velocity  components  v^  and  v^.  Note 

that  (4.14)  Is  linear,  and  that  (4.13)  Is  of  the  form: 

linear  terms  ■ ♦ non-1  inear  terms 

00 

so  that  as  M^-^0,  (4.13)  reduces  to  the  linear  (incompres- 
sible) continuity  equation. 

The  boundary  conditions  are  (2.3  - 2.10). 

If  the  direct  or  analysis  problem  is  solved,  2.10 


becomes 
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Vb  ®b 


v^b  sin  6^  = 0 


(4.15) 


where  and  Vj^  are  the  components  of  velocity  on  the 
body,  with  the  body  at  an  angle  9^^  measured  counter- 
clockwise from  the  z-axis. 

These  equations  are  said  to  be  in  primitive  variable 
form,  that  is,  they  are  still  in  terms  of  the  quan- 
tities of  interest,  the  velocity  components.  Note, 
however,  that  the  set  of  equations  (4.13)  and  (4.14) 
are  fully  equivalent  to  the  potential  equation  commonly 
used  for  inviscid  flows. 

The  pressure  coefficient  is  found  in  terms  of  the 
dimensionless  velocity  components  from: 


P-1 

I tmI 


where  P = [1 


t (Y-n 
2 


m!.(i  3 


y/(y-i  ) 


(4.16) 


In  the  limit  of  M^-^0  (incompressible  flow),  this 
becomes  : 

Cp  » 1 - (u"  + vM  (4.17) 

The  finite  el ement  formul ation  of  (4.13)  and  (4.14) 
requires  that  the  dependent  variables  v^  and  v^,  and 
for  axisymmetric  flow  the  independent  variable  r,  be 
approximated  across  each  element  by  interpolation 
functions  of  the  form 
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where  are  the  interpolation  functions, 

and  rj^  are  values  of  v^,  v^  and  r at  the  Nth  node, 

and  n is  the  number  of  nodes  per  element. 

The  simplest  possible  element  is  triangular, 
and  has  three  nodes  and  linear  interpolation  functions. 
Other  elements  may  have  other  shapes  and  higher  order 
interpolation  functions.  Each  element  in  the  flow 
field  gives  rise  to  a n*n  coefficient  matrix.  Thus, 
linear  triangular  elements  require  computation  of  a 
3*3  matrix.  Higher  order  interpolation  functions  may 
increase  accuracy,  but  also  require  lengthy  computa- 
tion of  large  coefficient  matrices.  Linear  triangular 
elements  assure  that  the  final  solution  will  have  con- 
tinuous values  of  the  dependent  variables  throughout 
the  flow  field  (called  C°  continuity).  Other  elements 
have  been  devised  that  also  assure  continuity  of  the 
first  m derivitives  of  the  variable  (called  C*"  con- 
tinuity). With  these  elements,  values  of  the  variable 
and  its  first  m derivitives  must  be  solved  for  at  the 
nodes,  increasing  solution  time. 

The  potential  equation  commonly  solved  by  finite 
difference  methods  is  second  order,  demanding  second 
order  elements  with  large  coefficient  matrices.  The 
solution  must  be  differentiated  to  obtain  the  velocity 
components,  making  continuity  desirable.  Both  re- 
quirements add  up  to  large  computer  times. 
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The  primitive  variable  formulation  given  by  equa- 
tions (4.13)  and  (4.14)  contains  only  first  derivitives, 
and  so  requires  only  first  order  elements  with  C°  con- 
tinuity. To  minimize  requirements  on  the  interpolation 
functions  and  thereby  reduce  computer  time,  the  primi- 
tive variable  formulation  of  subsonic  inviscid  flow  was 
chosen  over  the  potential  formulation. 

Isoparametric  quadrilateral  elements,  as  shown  in 
Fig.  4.5  employ  the  same  interpolation  functions  for  all 
variables  of  interest  in  arbitrarily  shaped  quadrilateral 
elements.  They  were  chosen  over  triangular  elements 
since  their  shape  is  more  suited  to  a roughly  rectangular 
flow  field,  and  since  their  interpolation  functions, 
while  almost  linear,  include  a cross  product  term  that 
increases  accuracy. 

Figure  4.5  shows  an  arbitrary  quadrilateral  element.  A 
non-dimensional  or  "natural"  coordinate  system  (e,n) 
is  established  at  the  centroid  of  the  element,  such  that 
the  coordinates  of  the  four  nodes  are  - 1.  This 
simplifies  the  resulting  expressions  for  the  interpola- 
tion functions.  Note  that  the  node  numbering  must  pro- 
ceed counterclockwise  around  each  element.  Otherwise, 
interpolation  functions  may  take  on  negative  values. 

It  is  impractical  to  derive  the  expressions  for 


the  interpolation  functions  and  their  derivitives  here. 
Details  are  available  in  [62]  or  [68],  and  all  results 
are  presented  here  for  convenience. 
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FIGURE  4.5  ISOPARAMETRIC  QUADRILATERAL  ELEMENTS. 

NON-DIMENSIONAL  (OR  "NATURAL")  COORDINATES 
(e  , n)  ARE  ESTABLISHED  AT  THE  CENTROID  OF 
THE  ELEMENT  SUCH  THAT  THE  COORDINATES  OF  THE 
NODES  ARE  t 1.  THIS  SIMPLIFIES  THE  EXPRES- 
SIONS FOR  THE  INTERPOLATION  FUNCTIONS.  NODES 
MUST  BE  NUMBERED  COUNTER  CLOCKWISE  AROUNG  THE 
1 ELEMENT. 
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Interpolation  for  any  Variable  <{>: 


0 = 


where 


Hiji  + + £!(})  + 514) 

l l 2 2 3 3 I*  4 


j (1-e)  (1-n) 
\ (l+e)(l-n) 
j (l+£)(l+n) 


I (1-£)(1+n) 


(4.1 9a) 


4 t J 

It  is  sometimes  convenient  to  rewrite  (4.19a)  in 
the  form: 

4)  s j [a  + be  + cn  + den] 
where 


a = 4>  ■^4>  "^4) 

12  3 4 

b»-(^  + <>  + (^  -1^ 

12  3 4 

c*-<^  - 4)  + 4)  + 4> 

12  3 4 

d“4^  ^4*  "^4^  *4^ 

12  3 4 


) (4.19b) 


Derivitives  of  Interpolation  Functions: 

The  finite  element  formulation  of  (6)  and  (7) 
requires  first  derivitives  of  the  interpolation  functions 


given  in 

(4.19a) 

. These  are  as  follows: 
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(4.20b) 
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The  Jacobian  matrix  [J]  relates  derivitives  in  the 
local  (C»n)  and  global  (z,r)  coordinate  systems. 


CJ] 
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(4.21a) 


and  its  determinant  |J|  is  given  by 


1 


where  are  functions  of  the  global  nodal  coordinates 


a » (z  -z  ) (r  -r  ) - (z  -z  ) (r  -r  ) 
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Integration 
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/ (4.21b) 


The  Galerkin  formulation  of  (4.13)  and  (4.14)  re 
quires  integration  of  various  functions  over  the  area 


no 


of  an  element.  The  relation  between  integrations  in 
global  and  local  coordinates  is: 

= /,]  /.|  | J | f (e,n)d£dn  (4.22a) 

Since  analytic  integration  is  often  impossible  due  to 
the  term  J , Gaussian  numerical  integration  is  used. 

L L 

/A,.oa^(z,r)d2dr  *11  u.-w.  | J | f ( e ,n ) (4.22b) 

Area  j k 

where  the  Gaussian  weight  functions  available 

elsewhere,  and  L is  the  order  of  the  Gaussian  integra- 
tion, (L=3  has  proven  to  be  sufficient). 

Galerkin  Formulation 

Now  the  approximations  for  v^,  v^  and  r given  by 
equation  (4.18)  are  substituted  into  the  equations  of 
motion  (4.13)  and  (4.14).  In  general,  the  right  hand 
sides  will  no  longer  equal  zero;  indeed,  they  will  each 
equal  some  residual.  The  Galerkin  method,  a subclass 
of  the  method  of  weighted  residuals,  is  used  to  minimize 
these  residuals. 

The  Galerkin  method  uses  the  interpolation  functions 
as  weight  functions,  and  requires  that  these  functions 
be  orthogonal  to  the  residuals  over  the  volume  of  any 
element.  Gartling  [75]  and  Oden  [76]  have  shown  that 
this  method  is  equivalent  to  an  integral  mechanical  energy 


balance;  so  the  method  has  physical  significance. 

Details  of  this  formulation  are  unnecessary.  The 
Galerkin  formulation  of  (4.13)  and  (4.14)  may  be  written 
in  the  following  matrix  form: 


C-M^[F(v  V ) + G(v,.v  )].(0+jE)  - M4CH(v  ,v  ) + jKv  ,v  )] 
z r z r zr  zr 

A , B 


(4.23) 

Here,  the  upper  row  is  the  Galerkin  formulation  for  the 
compressible  continuity  equation  (4.13)  and  the  bottom 
row  is  the  formulation  of  the  i rrotational i ty  condition 
(4.14). 

The  terms  A through  I are  each  (4*4)  coefficient 
matrices  given  by  the  following  integrals. 
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(4.24  continued) 
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2^  •^/\[(^l'^L^  ^ (^L^i^)  -l]Q|^f2nil  dzdr 


I = lllli 

^NM 


(4.24) 


M = 1 . 2 , 3 , 4 
M = 1 ,2,3,4 


This  set  of  algebraic  equations  is  non-linear  due 
to  the  terms  F,  G,  H,  and  I.  In  the  incompressible  limit 
of  = 0,  these  terms  vanish  and  the  system  becomes 
1 inear. 

A similar  set  of  equations  may  be  written  for  each 
element  in  the  flow  field.  Equations  from  each  element 
are  assembled  into  a global  matrix  equation  using  stand- 
ard techniques  available  in  [62]'  or  [68]. 

Dirichlet  boundary  conditions  (specified  values  of 
u or  V as  given  by  equations  2.7  and  2.8  can  now  be  sub- 
stituted directly  into  the  global  matrix.  Movement  of 
known  terms  to  the  right-hand  side  makes  the  right-hand 
side  non-zero  and  the  equation  set  non-singular. 

Neumann  boundary  conditions  (e.g.  the  tangency  con- 
dition given  by  equation  (4.15))  are  applied  via  LaGrange 
multipliers.  Details  may  be  found  in  [62]  or  [68]. 

The  tangency  constraint  is  applied  at  each  node  on  a 
body,  giving  rise  to  extra  equation  for  the  LaGrange 
multiplier  at  each  body  node.  Physically,  the  LaGrange 
multiplier  represents  the  "energy"  required  to  hold  the 
Neumann  boundary  constraint.  Practically,  the  value  of 
the  multiplier  is  useless. 


Iterative  Solution  of  the  Non-linear  Algebraic 
Equations 

Gartling's  work  with  a primitive  variable  form  of 
the  Navier-Stokes  equations  [75,77]  lead  to  the  follow- 
ing iterative  scheme  for  solving  the  non-linear  algebraic 
set  (4.23) 


(D+JE)  - Mf,[H(v""^  ,vjl‘^  " 

r n ") 

"z 

, •»/  n*l  n*l\T 
+ Jl(v^  ,v^  )] 

<t 

» 

B 

n 

- 

(4.25) 

where  the  superscript  ()^  refers  to  the  iteration  number. 

For  the  first  iteration  (n=l),  the  non-linear  terms  F,  G, 

H and  I are  set  to  zero,  and  the  resulting  linear  set  is 

solved  using  a Gauss-Jordan  scheme  for  banded  matrices. 

Thus,  the  n=l  solution  is  the  incompressible  solution, 

which  is  always  useful  for  comparison.  In  subsequent 

iterations,  the  non-linear  terms  are  calculated  from  the 

previous  (n-1)  values  of  u and  v. 

The  solution  is  considered  to  be  converged  when  all 

v"-v!]"^  and  all  v[I-v[I"^  are  less  than  0.0001.  This 
z z r r 

invariably  occurs  in  3 or  4 iterations,  unless  the  flow 
goes  locally  supersonic,  in  which  case  the  solution 
almost  never  converges. 
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4.3.3  Results 

In  each  of  the  following  examples,  an  automatic 
mesh  generation  subroutine  was  used  to  simplify  input 
to  the  program  and  to  minimize  the  band  with  of  the 
resulting  equations.  As  shown  in  Figure  4.6  this  sub- 
routine installs  nodes  along  vertical  columns  and  along 
roughly  parallel  rows.  If  the  number  of  nodes  per 
column  is  designated  Nr,  the  number  of  nodes  per  row 
is  Nz,  and  the  number  of  nodes  on  the  body  is  NB,  then 
the  total  number  of  nodes  is  NzxNr,  the  number  of  ele- 
ments is  ( Nz- 1 ) *( Nr- 1 ) , and  the  number  of  equations 
solved  is  2*Nz*Nr+NBo 

All  computer  times  given  below  are  for  an  IBM  370/153. 

Sphere 

Flow  over  a sphere  was  computed  on  a 9*27  node 
grid.  The  resulting  incompressible  pressure  coefficient 
is  shown  in  Figure  4.7  (circles),  compared  to  the  exact 
solution  (solid  line).  The  compressible  solution  at 
* 0.5  (plus  signs)  is  compared  to  the  Gothert's  rule 
compressibility  correction  (dashed  line)  on  the  same 
plot.  The  solutions  are  good  near  the  stagnation  region, 
but  worsen  near  the  peak  of  the  sphere  where  the  finite 
free  stream  boundary  has  the  most  effect.  Computer 
times  were  22  seconds  for  the  incompressible  solution, 
and  149  seconds  for  the  compressible  solution.  The 
compressible  solution  required  7 iterations  since  the 
flow  approached  the  critical  Mach  number  of  = 0.57. 


COUPLING  STUDY.  THERE  ARE  7*52  = 364  ELEMENTS,  (7+l)*(52+l)  = 424  NODES,  38  BODY  NODES 
AND  2*424+38  = 886  EQUATIONS.  THE  GRID  HAS  GENERATED  AUTOMATICALLY  BY  THE  COMPUTER 
PROGRAM. 
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FIGURE  4.7  and  4.8  COMPARISON  OF  FINITE  ELEMENT  AND  ANALYTIC  SOLUTIONS 

FOR  THE  PRESSURE  COEFFICIENT  ON  A SPHERE  AT  M^=0 
AND  M^=0.5  (4.7),  AND  ON  A CYLINDER  AT  M^=0  AND 
M = 0.38  (48). 
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The  arid  used  for  the  sphere  was  also  used  to 
compute  the  pressure  coefficient  over  a cylinder,  as 
shown  in  Figure  4.8.  Agreement  with  the  analytic  solu- 
tion is  similar  to  that  of  the  sphere,  but  slightly 
worse.  This  was  expected,  since  the  two-dimensional 
cylinder  represents  a larger  flow  disturbance  than  the 
axisymmetric  sphere.  Computer  times  were  22  seconds 
for  the  incompressible  solution  and  110  seconds  for  the 
M * 0.38  solution  (5  iterations). 

Rankine  Ovoid 

Figure  4.9  compares  the  calculated  and  exact  pres- 
sure coefficients  on  a 3.16:1  aspect  ratio  axisymmetric 

Rankine  ovoid.  At  M * 0,  the  finite  element  solution 

00 

(circles)  and  the  exact  solution  (solid  line)  agree 
almost  exactly.  At  = 0.6,  the  finite  element  solution 
(plus  signs)  compares  well  with  the  Gothert's  rule  compres- 
sibility correction.  A 10*29  node  mesh  was  used  over  the 
quarter  body.  Incompressible  and  compressible  solution 
times  were  29  seconds  and  120  seconds  respectively. 

Rankine  Oval 

Figure  4.10  compares  the  calculated  and  exact  pressure 
coefficients  on  a 3.16:1  aspect  ratio  two-dimensional 
Rankine  oval.  It  should  be  noted  that  the  body  profile 
is  not  the  same  as  that  of  the  ovoid  mentioned  above. 


FIGURE  4.9  and  4.10  COMPARISON  OF  FINITE  ELEMENT  AND  ANALYTIC  SOLUTIONS 


FOR  THE  PRESSURE  COEFFICIENT  ON  A 2-D  RANKINE  OVOID 

(4.9)  AND  ON  AN  AXISYMMETRIC  RANKINE  OVAL  (4.10), 

BOTH  AT  M = 0 AND  M = 0.6 
00  00 
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The  finite  element  incompressible  solution  (circles) 
and  the  exact  solution  (solid  line)  compare  well,  except 
in  the  region  of  peak  negative  pressure,  where  grid 
spacing  may  have  been  too  coarse.  A 9*29  node  grid 
was  used  over  the  quarter  body.  The  = 0.6  solution 
is  also  shown  (plus  signs).  Computer  times  were  25 
seconds  at  = 0 and  123  seconds  at  M = 0.6. 

14%  Joukowsky  Airfoil 

The  finite  element  program  has  been  used  to  solve 
for  the  flow  over  a 14%  thick  symmetric  Joukowsky  airfoil 
at  = 0 and  = 0.6.  Figure  4.11  shows  the  exact 
pressure  coefficient  (solid  line),  the  incompressible 
finite  element  solution  (circles),  the  = 0.6  finite 
element  solution  (plus  signs),  and  the  M = 0.6  Gothert's 
rule  compressibility  correction  to  the  exact  solution. 
Both  the  incompressible  and  compressible  solutions  show 
excellent  agreement  with  the  analytic  solutions.  The 
only  error  is  again  in  the  peak  negative  pressure  region, 
probably  indicating  inadequate  mesh  spacing.  The  mesh 
consisted  of  7 node  rows  by  37  node  columns,  for  a total 
of  259  nodes  (13  on  the  body).  This  amounts  to  216 
elements  and  531  simultaneous  equations.  The  incompres- 
sible solution  took  21  seconds,  and  the  M * 0.6  solu- 

00 

tion  took  83  seconds  (4  iterations).  These  times  are 
not  out  of  line  with  finite  difference  methods. 
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Axlsymmetric  Boattail  Model 

Chow,  Bober  and  Anderson  [60]  published  experi- 
mental data  and  finite  difference  calculations  for  the 
pressure  coefficient  on  a NASA  axisymmetric  boattail 
model.  Figure  4.5  shows  the  B^BS  node  finite  element 
grid  used  to  re-compute  this  flow  in  the  present  study. 

Figure  4.12  shows  the  incompressible  finite  element 
solution  (solid  line)  compared  to  a surface  source  method 
(dashed  line).  Qualitative  agreement  is  good,  but 
quantatively  there  is  some  disagreement  in  the  results. 

It  is  thought  that  this  is  due  to  the  finite  location 
of  the  free-stream  boundary.  This  boundary  is  30  body 
radii  away,  but  only  1.6  body  lengths  (excluding  the 
sting)  away.  A compressible  = 0.8  solution  is  also 
presented  (circles). 

Boundary  Layer  Coupling 

Sasman  and  Cresci's  compressible  turbulent  boundary 
layer  program  [61]  was  coupled  iteratively  to  the  finite 
element  inviscid  flow  prgoram  using  the  classical  method 
of  augmenting  the  body  by  the  displacement  thickness. 

Briefly,  Sasman  and  Cresci  reduce  the  integral 
momentum  and  moment  of  momentum  equations  to  an  "incompres- 
sible" form  via  a Mager  type  transformation.  They  use 
power  law  velocity  profiles  and  the  Ludwei g-Ti 1 1 man 
skin  friction  correl 1 ation . Their  dependent  variables 


FIGURE  4.12  COMPARISON  OF  FINITE  ELEMENT  AND  SURFACE  SOURCE  SOLUTIONS  FOR  THE 
PRESSURE  COEFFICIENT  ON  A NASA  BOATTAIL  MODEL,  AND  THE  COMPRESSIBLE 
FINITE  ELEMENT  SOLUTION  FOR  M =0.8. 
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are  the  momentum  thickness  9 and  the  shape  factor 
H = 6*/e,  where  6*  Is  the  displacement  thickness. 

Chow,  Bober  and  Anderson  [60]  used  this  method  for  their 
boattail  work  at  Mach  numbers  from  .56  to  .9,  with  good 
results  at  all  but  the  highest  Mach  number. 

Figure  4.13  compares  the  pressure  coefficients  over 
the  boattail  model  at  M = 0.8  calculated  by  the  inviscid 
method  alone  (circles),  and  by  the  i nv i sci d/ vi scous  inter- 
action (solid  line). 

Figure  4.14  is  an  enl argement  of  the  boattail  region 
from  Figure  4.13  compared  with  the  results  of  Chow,  Bober 
and  Anderson  [60].  The  finite  el ement/Sasman-Cresci 
method  underestimates  the  pressure  coefficient  in  this 
region.  It  should  be  noted,  however,  that  the  finite 
element  grid  for  this  problem  was  8 nodes  high  by  53  nodes 
long,  while  the  grid  used  by  Chow,  et  al . was  26  nodes 
high  and  101  nodes  long.  Furthermore,  Chow,  Bober  and 
Anderson  used  an  infinite- to -finite  transformation  on 
the  flow  field  in  the  r direction,  while  the  present  study 
used  a finite  free-stream  boundary.  Chow,  Bober  and 
Anderson's  finite-difference  interaction  problem  took  7 
minutes  on  a CDC  6600,  while  the  finite  element  method 
took  13  minutes  on  an  IBM  370/158. 

Finally,  Figure  4.15  compares  the  displacement  thick- 


i 


nesses  in  the  boattail  region  calculated  by  Chow,  et  al  . 
and  by  the  present  study.  The  disagreement  is  undoubtedly 
due  to  discrepancies  in  the  inviscid  solutions. 


FIGURE  4.13  COMPARISON  OF  PRESSURE  COEFFICIENTS  ON  A NASA  BCATTAIL  MODEL  AT 


M^=0.8.  THE  "+"  SIGNS  ARE  THE  INVISCID  SOLUTION  ON  THE  ACTUAL 
BODY.  THE  SOLID  LINE  IS  THE  INVISCID  SOLUTION  ON  THE  BODY  AUG- 
MENTED BY  THE  DISPLACEMENT  THICKNESS  6*  CALCULATED  ITERATIVELY 
FROM  A SASMAN-CRESCI  TURBULENT  BOUNDARY  LAYER  DECK  AT  M„=  0.8, 

R^- 3.59*10'^ 
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FIGURE  4.U  COMPARISON  OF  FINITE  ELEMENT  AND  FINITE  DIFFERENCE 
SOLUTIONS  FOR  THE  PRESSURE  COEFFICIENT  ON  A NASA 
BOAHAIL  MODEL  AT  M„=«0.3.  BOTH  ARE  SHOWN  WITH  AND 
WITHOUT  BOUNDARY  LAYER  INFLUENCE,  CALCULATED  USING  A 
SASMAN-CRESCI  TURBULENT  BOUNDARY  LAYER  DECK. 
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4.3.4  Design  or  Inverse  Problem 


Inverse  boundary  layer  coupling  requires  that  the 
inviscid  flow  method  be  capable  of  determining  the  edge 
angle  9 given  the  edge  Mach  number  or  the  total  edge 
velocity  V^.  This  problem  is  roughly  equivalent  to  the 
classic  design  problem  in  which  a body  geometry  is  to  be 
determined  from  a specified  pressure  coefficient  C^. 

Design  calculations  were  attempted  with  the  finite 
element  program  using  the  following  iterative  scheme. 

An  initial  guess  is  made  at  the  body  geometry,  and 
hence  the  edge  angle  9.  Assuming  that  the  flow  follows 
the  initial  geometry,  the  pressure  coefficient,  edge 
Mach  number,  or  total  edge  velocity  may  be  algebraically 
resolved  into  v^  and  v^  velocity  components  on  that  body. 

In  an  analysis  problem,  the  Neumann  boundary  condi- 
tion of  flow  tangency  is  applied  on  the  body.  In  the 
inverse  or  design  problem,  that  condition  is  replaced  by 
a Dirichlet  boundary  condition  in  which  one  of  the  two 
velocity  components,  either  v^  or  v^,  is  specified  on  the 
initial  geometry.  The  finite  element  method  solves  for 
the  other  component. 

Since  the  geometry  guessed  initially  probably  does 
not  correspond  to  the  desired  distribution  of  Cp,  M^,  or 
V^,  the  velocity  component  specified  and  the  velocity 
component  solved  for  will  not  satisfy  the  requirement  of 
flow  tangency.  Physically,  the  solution  represents  a 


body  with  bl owi ng  or  sucti on , which  displaces  the  stagnation 
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streamline  from  the  guessed  body. 

The  displacement  of  the  stagnation  streamline  may 
be  used  as  a basis  for  updating  the  guessed  geometry. 

L.  A.  Carlson  [47]  developed  a finite-difference  design 
method.  In  his  work,  he  found  the  location  of  the  stag- 
nation streamline  by  integrating  the  expression 


In  =»  In 

dz  ~ v_ 


Guessed  Body 


then  used  the  location  of  this  streamline  as  an  updated 
body  geometry  for  the  next  iteration. 

T.  L.  Tranen  [46]  developed  a design  method  that 
used  a more  complex  update  technique.  This  technique  was 
based  on  adding  to  the  guessed  body  a kind  of  displace- 
ment thickness  of  sufficient  height  to  exactly  carry 
the  normal  mass  flow  resulting  from  the  blowing  or 
section.  During  the  course  of  the  present  study,  it  was 
discovered  that  certain  terms  in  Tranen's  scheme  were 
insignificant,  so  that  his  update  scheme  reduced  numer- 
ically to  Carlson's. 

One  advantage  of  the  finite  element  method  over 
finite  difference  methods  in  this  problem  is  that  inter- 
polation functions  are  available  for  the  terms  v^  and 
v^.  Hence,  the  differential  equation  for  the  updated 
geometry  can  be  evaluated  analytically  from  node  to 
node . 


Initially,  a small  region  of  the  boattall  model 


was  to  be  "designed"  from  a specified  pressure  dis- 


tribution determined  from  a previous  finite  element 


analysis.  To  test  the  method,  the  exact  corresponding 


geometry  was  Input  as  the  first  guess. 


When  the  velocity  components  were  applied  as 
the  Dirichlet  boundary  conditions  on  the  region  to  be 


"designed",  the  finite  element  method  behaved  unstably 


and  returned  a physically  unrealistic  flow  solution. 


No  explanation  Is  currently  available  for  this  phenomena. 


except  that  the  method  Is  numerically  unable  to  handle 


this  particular  boundary  condition. 


When  the  v^  velocity  components  were  specified,  the 


method  returned  the  correct  solution  for  the  v^  compon- 


ents. Thus,  this  approach  appears  promising.  When  the 


Initial  body  guess  was  not  the  body  corresponding  exactly 


to  the  specified  pressure,  the  method  still  returned 


reasonable  values  of  v^  In  the  region  to  be  designed. 


However,  the  update  scheme  tended  to  diverge  from  the 


true  body< 


Conci usions 


Experience  with  the  finite  element  formulation  of 


Inviscid  compressible  flow  In  primitive  variables  Indi- 


cates that  the  analysis  method  Is  both  fast  and  accurate 


The  ability  of  the  method  to  conform  to  arbitrary  bound- 


method'  in  the  analysis  problem  and  particularly  in 
the  iterative  problem  of  boundary  layer  coupling.  For 
the  Inverse  or  design  problem,  the  method  appears  promis- 
ing, but  further  investigation  is  needed. 

Recommendations  for  Further  Research 

Further  research  should  be  conducted  on  use  of  the 
finite  element  method  as  an  inviscid  flow  solution  tech- 
nique. The  method  should  be  tested  on  finer  grids  and 
faster  machines. 

The  problem  of  a finite  free-stream  boundary  could 
be  overcome  using  an  i nf i n i te- to-f 1 n i te  transformation 
of  the  r-coordi nate.  This  approach  was  used  by  Chow, 
et  al.,  in  their  finite  difference  boattail  study  [60]. 
Habashi  developed  a semi - i nf i n i te  boundary  elonent  [74] 
that  appears  to  be  more  consistent  with  the  present  work. 

Convergence  problems  of  the  design  or  inverse  method 
need  to  be  explored  and  eliminated,  so  that  inverse 
boundary  layer  coupling  may  be  tested  in  fully  compres- 
sible flow. 

Finally,  extensions  of  the  analysis  method  to  lift- 
ing surfaces  and  three-dimensional  bodies  would  be  highly 
des i rabl e . 


5. 
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Combination  of  the  Boundary  Layer  and  Inviscid  Flow 
Methods  Into  a Complete  Separated  Flow  Prediction 
Method 

Of  the  inviscid  flow  methods  discussed  in  Chapter  4, 
only  the  surface  source  method  has  thus  far  been  developed 
to  the  point  of  being  ready  to  combine  with  the  (successful) 
boundary  layer  method  of  Chapter  3.  In  this  chapter,  exper- 
iences with  the  attempt  to  combine  the  two  methods  will  be 
reported . 


5.1  Flow  Geometry  Selected 

The  flow  geometry  which  was  selected  for  the  at- 
tempts at  a complete  calculation  of  a flow  exhibiting 
separation  is  shown  in  Figure  4.3.  Separation  is  ex- 
pected to  take  place  on  the  boattail,  with  subsequent 
reattachment  on  the  sting.  This  geometry  is  typical 
of  that  employed  to  study  nozzle  afterbody  drag  and 
considerable  effort  has  been  expended  on  research  into 
the  development  of  an  adequate  method  for  calculating 
the  flow  over  such  a body  [6,45,60]  with  a view  to 
analytical  predictions  of  afterbody  drag.  The  separa- 
tion bubble  occuring  in  such  a flow  is  relatively  thin 
and  short  and  may  be  expected  to  be  steady.  It  was 
assumed  that  the  surface  was  smooth.  Air  at  a free 
stream  reference  Reynolds  number  (c„/v„)  of  1.7  x 10®(-r-^ — r'; 
was  assumed  to  approach  the  body  at  Mach  numbers  of 
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0.5  to  0.8.  It  was  expected  that  separation  would 
occur  earlier  at  higher  Mach  numbers;  this  was  verified 
by  the  calculations. 

5.2  Geometric  Problems  Involved  in  Combining  the 

Boundary  Layer  and  Inviscid  Flow  Methods 

Several  details  must  be  worked  out  in  order  to 
successfully  combine  the  methods  and  mak3  iterative 
calculations.  The  source  of  most  of  the  trouble  is 
due  to  three  items: 

(1)  The  boundary  layer  method  is  formulated 
in  a coordinate  system  parallel  to  and 
perpendicular  to  the  body  surface,  while 
the  inviscid  flow  procedure  is  formulated 
in  a coordinate  system  with  axes  along  an 
axis  of  symmetry  and  perpendicular  to  it. 

(2)  The  boundary  layer  thickness  is  added  on 
perpendicul ar  to  the  body  to  define  the  "6 
surface".  The  "edge  Mach  number"  and 
angle  Oat  a surface  point  x are  actually 
the  values  above  z(x)  + Az  a distance  of 

6 cos  a (see  Figure  5.1). 


* The  assumption  that  dominant  cross  stream  coordifi^te 
"y"  is  perpendicular  to  the  surface  is  of  less  validity 
in  the  neighborhood  of  separation,  however,  it  is 
nevertheless  retained. 


'Az  ' 

FIGURE  5.1  LOCATION  OF  POINT  AT  EDGE  OF  BOUNDARY  LAYER  IN  z,  r 
COORDINATE  SYSTEM 
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(3)  The  adding  of  a finite  6 perpendicular  to 
the  surface  results  in  an  ''overlapping"  of 
the  6 surface  in  the  immediate  vicinity  of 
a concave  corner  such  as  the  boattai 1 -sti ng 
j uncture . 

The  latter  two  of  these  items  are  a consequence 
of  the  first.  Item  (2)  is  simply  handled  by  the  trans- 
formation 

Mg(x,6)  Mg(z+A2,R+Ar ) 

0 (x,6)  + a{  x)-»-^  9 ( z+Az  , R+Ar ) 

^ _ J ^ i 

"Bounday  Layer  "In viscid  Flow 
Coordinates"  Coordinates" 

where 

Az  » 6{x)  sin  o(x) 

Ar  * 5 ( X ) cos  ot( X ) 

5 from  the  most  recent  boundary  layer  calculation  is  used. 

Item  (3)  is  more  complicated  since  it  involves  a 
physically  impossible  situation.  Nakayma  et  a1.  [9] 
introduced  a method  of  overcoming  a similar  difficulty 
by  defining  a triangular  control  volume  with  the  concave 
corner  at  one  vertex  and  with  sides  perpendicular  to  the 
upstream  and  downstream  surfaces.  The  boundary  layer 
equations  are  solved  up  to  the  inlet  face  of  the  control 
volume  and  terminated.  The  calculations  are  restarted 
at  the  downstream  face  with  parameters  determined  by  rrass 
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and  momentum  balances  for  the  control  volume.  A 
similar  method  was  investigated  for  the  present 
purposes,  employing  a control  value  composed  of  two 
congruent  triangles  as  shown  in  Figure  5.2.  This 
approach  was  finally  considered  too  complicated  and 
time  consuming  for  the  accuracy  gained.  The  treatment 
finally  employed  simply  marches  the  boundary  layer 
calculations  through  this  region  (the  "corner"  is  not 
recognized  in  the  surface  oriented  boundary  layer  coor- 
dinate system).  Points  in  the  region  of  interference 
(that  is  points  between  the  tail  location,  z.^.,  and 
Zj + sin  dj)  are  simply  discarded  when  transforming 
to  the  inviscid  flow  geometry.  In  this  way,  a smooth, 
continuous  "5  surface"  was  defined. 

5.3  The  Iterative  Procedure  and  Initial  Guesses 

Because  the  boundary  layer  equations  are  parabolic 

★ 

and  the  inviscid  flow  equations  are  elliptic  , simul- 
taneous solution  of  the  boundary  layer  and  inviscid 
flow  equations  is  not  possible  (as  it  is  in  supersonic 
separated  flow  calculations  [35-37])  and  iteration 
between  a complete  boundary  layer  calculation  and  a 
complete  inviscid  flow  calculation  must  be  employed. 

The  critical  variables  of  the  iteration  are  the  fluid 


* Thus  preserving  the  well  known  elliptic  nature  of 
the  separated  flow  problem. 
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velocity  (or  Mach  number)  and  angle  along  the  boundary 
between  the  viscous  ("boundary  layer")  region  and  the 
external  inviscid  flow  region.  The  calculation  pro- 
cedure in  each  region  calculates  "new"  values  for  one 
of  these  variables  at  each  point,  which  is  handed  to 
the  procedure  for  the  other  region  where  the  opposite 
variable  Is  updated  and  so  on  until  convergence  is 
achieved. 

There  are  several  possibilities  for  starting  the 
calculations.  The  most  obvious  is  to  first  calculate 
the  inviscid  flow  over  the  body  assuming  there  is  no 
boundary  layer  (a  "direct"  or  analysis  calculation)  and 
then  use  the  pressure  (velocity)  distribution  thus 
calculated  as  input  to  the  boundary  layer  method  to 
update  the  geometry,  etc.  Another  obvious  alternative 
is  to  first  calculate  the  boundary  layer  on  the  body 
assuming  no  "effect"  of  the  inviscid  flow,  that  is 
calculate  the  boundary  layer  under  the  assumption  of 
constant  pressure  and  then  use  the  calculated  edge 
angle  and  "5-surface"  as  input  to  the  inviscid  flow 
method  and  so  on.  In  actual  practice  these  approaches 
are  essentially  identical;  the  "constant  pressure" 
boundary  layer  does  not  significantly  modify  the  body 
shape  so  that  the  "bare  body"  pressure  distribution 
results  from  the  inviscid  calculation  in  the  second 
case.  Thus  in  either  case,  the  calculations  "start" 
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with  a pressure  distribution  essentially  identical  to 
the  value  for  the  completely  inviscid  flow,  and  the 
first  significant  problems  develop  in  the  boundary 
layer  routine. 

The  boundary  layer  calculations  are  started  in  the 
weak  interaction  mode  and  continue  in  this  mode  up  to 
a position  about  midway  down  the  boattail,  where  a 
(typical)  sharp  rise  in  6,  <5*,  H,  and  ^ occur.  These 
are  of  course  an  indication  of  impending  boundary 
layer  separation.  As  soon  as  one  of  the  criteria  of 
equations  3.34  are  satisfied  (typically  both  are  satis- 
fied at  the  same  calculation  step),  the  calculations 
switch  to  the  strong  interaction  mode.  A major  problem 
is  immediately  encountered  in  that  the  strong  inter- 
action mode,  0 is  required  as  input,  with  to  be  cal- 
culated; thus  a 0 distribution  must  be  assumed  to  con- 
tinue the  calculations.  In  order  to  obtain  a realistic 
boundary  layer  calculation,  the  O di stribution  must  be 
realistic.  The  distribution  must  have  the  general 
shape  shown  in  Figure  5.3.  The  rapid  Increase  of© 
associated  with  separation  must  be  maintained  for  a dis- 
tance but  the  value  of  © must  not  be  too  large  or  numer- 
ical problems  will  occur.  In  order  to  "turn"  the  flow 
back  toward  the  sting  and  "force"  a reattachment  €> 
must  decrease  and  attain  negative  values.  As  the  bound- 
ary layer  is  expected  to  redevelop  in  typical  fashion 
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after  reattachment,  an  eventual  return  to  positive  (small) 
O values  is  anticipated.  After  a considerable  amount 
of  "cut  and  try",  the  following,  procedure  for  con- 
structing a trial  Q distribution  was  settled  on.  Presz 
and  Pitkin  [6]  have  shown  that,  for  flows  separating 
from  boattailed  afterbodies  the  reattachment  point  can 
be  approximately  located  by  extending  a straight  dividing 
streamline  from  the  separation  point  to  its  intersection 
with  the  sting.  The  angle  this  line  makes  with  the  body 
surface  at  separation  is  a function  of  the  local  Mach 
number.  Using  this  method,  an  approximate  reattachment 
point,  Zj^,  is  determined. 

In  order  to  obtain  a realistic  (increasing)  Q 
distribution  downstream  of  separation,  it  is  assumed 
that,  up  to  the  boattail  -sting  junction,  the  velocity 
vector  is  parallel  to  the  axis  of  symmetry,  thus  0 is 
the  negative  of  the  body  slope  and  reaches  a maximum 
(in  this  case  about  .6  radians)  above  the  boattail  -sting 
junction.  Downstream  of  this  point,  Q varies  according 
to 

R T 

Thus  Q is  a minimum  at  the  estimated  reattachment  point. 
When  Q reaches  a value  of  zero  downstream  of  Zj^,  Q is 


continued  at  zero. 
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This  distribution  of  O is  sufficient  to  continue 
calculations  downstream  of  the  point  where  the  switch 
from  weak  to  strong  interaction  occurs;  however,  the 
values  of  O calculated  by  the  weak  interaction  formula- 
tion at  the  last  few  stations  before  the  switch  are 
usually  too  large;  often  larger  than  O max.  Those  values 
of  ©which  are  greater  than  0.3  radians  are  replaced  by 
values  computed  from  a second  order  curve  fit  between 
©=  0.3  and  the  initial  value  of  0(=-a)  at  the  switch 
point. 

Having  completely  specified  O,  calculations  are 
continued  in  the  strong  interaction  mode;  typical  re- 
sults are  shown  in  Figure  5.4.  At  some  point  fol- 
lowing reattachment,  the  conditions  of  equations  3.35 
are  satisfied  and  the  calculations  are  switched  back  to 
the  weak  interaction  mode  (M^  specified).  The  values  of 

needed  to  continue  the  calculation  are  available  from 
e 

the  previous  (initial)  inviscid  calculation;  however, 
it  usually  occurs  that  the  last  computed  from  the 
strong  interaction  mode  is  not  equal  to  the  inviscid 
value  at  that  point  so  a discontinuity  results.  This 
is  removed  by  an  exponential  fairing  between  the  last 
strong  interaction  value  and  the  value  of  five 
boundary  layer  thicknesses  downstream.  Calculations 
are  then  continued  to  the  tail  of  the  sting  far  downstream. 


The  results  of  the  boundary  layer  calculation  are 
a set  of  values  of  6(x),  "new"  values  of  0(x)  at  points 
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where  weak  interaction  calculations  were  made,  and 
"new"  values  of  Mg(x)  at  points  where  strong  inter- 
action calculations  were  made.  The  inviscid  flow  pro- 
cedure is  now  set  up  to  find  "new"  M^'s  for  the  weak 
interaction  points  (from  the  newly  cal cul ated  G' s ) and 
"new"  O' s for  the  strong  interaction  points.  The  "6 
surface"  is  constructed  and  the  equations  for  the  re- 
quired Mg  and  0 (actually  V/V^  and  9)  values  formulated, 
as  outlined  in  Section  4.2.  As  an  initial  guess  to 
start  the  Newton  iteration  of  equations  4.6  - 4.8,  all 
source  strengths  are  set  equal  to  zero,  all  unknown 
edge  velocities  are  set  equal  to  the  free  stream  value 
(V/V^  = 1)  and  unknown  angles  (0)  are  set  equal  to  zero. 

At  this  point,  the  method  breaks  down.  All  complete 
calculations  attempted  to  date  have  failed  to  converge 
in  the  Newton  iteration  in  the  inviscid  flow  program. 
Characteristically,  the  failure  to  converge  is  dominated 
by  the  following  problems: 

The  values  of  sin  9 computed  at  the  "strong  inter- 
action" points  tend  to  oscillate  between  positive  and 
negative  values  from  point  to  point  and  iterations  to 
iterations.  If  not  strongly  controlled,  the  values  of 
sin  9 will  exceed  1.0.  (See  Figure  4.4.) 

The  cycle  to  cycle  error  (monitored  by  the  average 
change  of  either  V/V  ).  or  sin  o).  between  iterations 
decreases  to  a minimum  of  about  4 iterations  then  in- 
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creases.  The  relative  minimum  error  is  rather  large. 

The  "solutions"  in  the  neighborhood  of  switch  points 
are  discontinuous. 

Several  "obvious"  remedies  to  alleviate  this  be- 
havior were  tried;  these  include: 

(1)  Undercorrection  using  only  a portion  of  the 
newly  calculated  (from  the  boundary  layer 
method)  5,  O combined  with  the  previous 
values  e.g. 

*^new  * f *^old  I ^new 

(2)  "Smoothing"  of  the  6 surface  and  the  Q 
distributions. 

(3)  Using  the  flow  over  the  "bare  body"  or  the 
(assumed  solid)  "5  surface"  as  the  initial 
guess  for  the  Newton  iteration. 

None  of  these  "fixes"  was  observed  to  improve  the  situ- 
ation. 

A possible  explanation  and  resolution  of  the  prob- 
lem is  as  follows.  By  considering  the  boundary  layer 
calculation  method,  one  may  conclude  that  although  calcu- 
lations are  possible,  one  must  tread  a very  tortuous 
path  to  make  them.  In  developing  the  boundary  layer 
method,  it  was  assumed  that  whatever  was  needed  to  con- 
tinue could  be  obtained  i.e.  mode  switch  was  possible 
where  necessary,  singularities  were  carefully  avoided,  etc. 
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Using  measured  data  with  the  boundary  layer  method 
verified  its  ability  to  make  calculations  if  "good" 
information  is  available  as  input.  In  a complete  calcu- 
lation, since  all  liberties  have  been  taken  with  the 
boundary  layer  method,  the  i'nviscid  flow  method  must 
take  up  the  slack,  making  exactly  the  calculations  re- 
quired by  the  boundary  layer  method.  We  cannot  "play" 
with  the  inviscid  routine  to  find  out  what  types  of 
boundary  conditions  are  most  efficient  or  acceptable  at 
various  points. 

It  is  still  felt  that  the  method  is  capable  of 
making  complete  separated  flow  calculations  but  that 
the  inviscid  flow  method  will  converge  only  if  called 
upon  to  "improve"  or  Q values  which  are  close  to 
actual  values;  in  other  words,  an  initial  guess  close 
to  the  final  result  must  be  provided.  Two  methods  are 
suggested  for  providing  this  initial  guess.  Firstly,  a 
semi -empi rical  approximate  method  of  estimating  the  flow 
parameters  (especially  pressure  distribution)  such  as 
those  proposed  by  Presz  and  Pitkin  [6]  or  Kuhn  [78] 
could  be  used  to  generate  an  initial  guess  for  the  pres- 
ent method. 

Secondly,  an  interactive  computer  terminal,  hope- 
fully with  graphics  capability,  could  be  employed.  The 
"man  in  the  loop"  could  postrate  various  reasonable 
initial  guesses  until  one  which  allowed  calculations  to 
proceed  is  found. 
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At  the  time  this  research  was  carried  out,  the 
latter  possibility  was  not  available  to  the  authors. 
Toward  the  end  of  the  project,  attempts  were  made  to 
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persue  a method  of  the  first  type,  similar  to  the  method 
of  Kuhn,  in  which  a polynominal  distribution  of  6 is 
assumed  in  the  strong  interaction  region,  with  coef- 
ficients selected  in  order  to  minimize  the  deviation 
of  the  calculated  edge  velocities  from  the  prescribed 
values.  No  formal  procedure  for  constructing  the  poly- 
nominal  approximation  to  minimize  the  deviation  was 
programmed;  efforts  to  select  a set  of  coefficients  by 
"cut  and  try"  were  finally  abandoned  as  too  time  con- 
sumi ng . 

At  the  present  time,  it  is  hoped  that  the  finite 
element  procedure,  still  under  development,  will  be  capable 
of  combination  with  the  boundary  layer  method  for  a com- 
plete calculations  routine. 
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6.  Conclusions  and  Recommendations  for  Further  Research 

The  following  are  considered  to  be  the  major  accomplish- 
ments of  this  work  and  conclusions  that  may  be  drawn  from  it. 

(1)  The  strong  interaction  formulation  is  a viable 
framework  for  separating  boundary  layer  analysis. 

(2)  The  "boundary  layer  equations"  are  adequate 
representations  of  the  flow  in  the  viscous  region. 

(3)  It  is  possible  to  develop  a general  integral 
method  for  solving  the  boundary  layer  equations 
in  cases  involving  flow  separation. 

(4)  Simple  algebraic  (non-equilibrium)  turbulence 
models  are  adequate  for  use  with  such  a method 
and  the  effects  of  "normal  stresses",  both  in 
the  mean  flow  equations  and  in  the  turbulence 
model  are  not  significant. 

(5)  The  integral  boundary  layer  method  is  especially 
prone  to  the  appearance  of  "singularities". 

These  singularities  are  usually  mathematical 
rather  than  physical  and  smooth  passage  through 
them  should  not  be  a forced  condition  on  the 
analysis.  The  singularity  can  usually  be  avoided 
by  switching  between  the  "strong"  and  "weak" 
interaction  modes  of  calculation. 
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(6)  The  integral  method  probably  does  not  offer  any 
advantage  over  a finite  difference  or  finite 

4 

element  formulation  in  terms  of  ease  of  formula- 
• tion,  ease  of  programming,  and  computation  time.  | 

The  advantages  of  the  method  thus  lie  in  the  fact 
that  the  desired  Interaction  parameters  (M^, 

©,  5)  are  explicit  in  the  formulation,  and  it 
seems  rather  forgiving  in  terms  of  the  turbulence 
model  employed. 

(7)  Experience  with  the  method  of  integral  relations 
for  the  inviscid  flow  indicates  that  the  method 
is  too  complex  and  several  of  the  required  steps 
too  artificial  for  the  method  to  be  given  serious 
consideration. 

(8)  It  is  possible  to  develop  an  a priori  "inverse" 
or  "design"  calculation  procedure  for  inviscid 
flow  over  axisymmetric  bodies  using  the  surface 
singularity  method. 
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(9)  The  ability  to  make  complete  a priori  separated 
flow  calculations  is  apparently  restricted  by  the 
need  for  initial  guesses  close  to  the  final  answer. 

(10)  The  finite  element  method  is  a viable  alternative 
to  the  finite  difference  method  for  computing 
compressible  (non-linear)  flows  in  the  "direct" 
or  "analysis"  mode  and  can  be  successfully  com- 
bined with  boundary  layer  calculations  via 
"displacement  thickness"  interaction. 


i 


1 50 

(11)  Problems  still  remain  in  the  development  of  the 
finite  element  approach  to  the  "inverse"  or 
"design"  problem  of  inviscid  flow. 

The  following  are  recommendations  for  further  work 
suggested  by  the  results  of  this  study. 

(1)  The  development  of  an  "inverse"  finite  element 
procedure  should  continue. 

(2)  A procedure  such  as  that  of  Presz  and  Pitkin  [6] 

Kuhn  [78]  should  be  incorporated  to  obtain  an 
initial  estimate  of  the  flow  parameters,  which 
can  be  improved  by  the  present  method. 

(3)  Apparently,  the  use  of  curved  rather  than  straight 
body  elements,  polynominal  rather  than  uniform 
singularity  distributions,  and  combined  singular- 
ities rather  than  simple  sources  offers  consider- 
able improvement  in  the  speed,  stability,  and 
accuracy  of  the  surface  source  method,  especially 
for  the  design  problem  [56-58].  Accordingly, 
such  improvements  should  be  introduced  into  the 
current  method. 

(4)  "Displacement  thickness"  interaction,  especially 
in  the  strong  interaction  calculation  ((S*  rather 
than  © specified  in  the  boundary  layer  equations) 
should  be  investigated  for  the  present  program. 
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(5)  A finite  element  formulation  of  the  boundary 
layer  flow  should  be  Investigated.  In  this 

i 

(differential)  framework,  the  algebraic  turbulence 
model  should  be  compared  with  one  and  two  equa- 
tion models. 

(6)  Finite  difference  models  of  both  the  boundary 
layer  and  inviscid  flow,  retaining  the  strong 
interaction  formulation  of  Chapter  2,  should  be 
investigated. 
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Below,  the  various  derivatives  of  the  velocity  profll 
equations  3.5  and  3.7  are  tabulated. 


3a 


(y-I)RF 


" yVMjl-9/4 


(y-DM^  RF(Y-l)Mi 

, v.i°  , - 7/4  Trf-^] 


■ 6*/M  [1-9/4  - 7/4 

='"e  ® 1 +1^ 


(y-DM? 


RF(y-i )m^ 

1 +rf(1^)m; 


If  - /TTV  c|  In  d+y'")  + B - (l.Sy'’  + B )e' * 

If  * /T^Up'  [fin  (1+5'^)  + B - (1  + B)e"'^®‘^^] 


1^  “ y*/^ 


IL.  = A*/T, 


M , . /TTv  [f  C Ug  sin  (|c  ) cos  (f  E )] 


+ (1  .55*  + B - 8.333)] 


/mV  [ 

k(l+5*) 


/TTV  [XC W-  + .18e"'''®y  (1.5y*  + B-8 

k(l+y  ) 


'l-aH*  CX6^( ^r-  + .18e'*^®^  {1.5y*  + B 

kd+y*) 


► iTUg  sin  (jC  ) cos  (jC  )] 


/nTv  sin*  (5  C ) 


JT- 


2a2 


a‘((» 


(y-I)RF  m"  D r* 


- <|)/a  + /--  ( 

7 a 


s in~^  a<(i^ 


. F/a  ♦ yTITT*  In  (1+5*)  + B - ( 1 . 55*  + 8 )e' 


i^M2^-7/4 


R5MJXI  (1  +1^  +RF(^)M^) 


.333)] 

3.333)) 


.185*  ^ 


155 

List  of  Symbols 

a Mach  number  function  defined  in  Appendix  A 
8 "Law  of  the  wall"  constant 

c Speed  of  sound 

Cx  Skin  friction  coefficient  C-  = p^uf 

T f w £ e e 

F Ratio  of  shear  to  normal  stress  turbulence  production 

H Boundary  layer  form  factor  H ■ 6*/9 

★ 

Boundary  layer  form  factor 

j * 0 for  plane  2-D  flow,  » 1 for  axisymmetric  flow 
|J|  Jacobian  determinant 

l Mixing  length 

M Mach  number 

P Pressure 

q Turbulence  kinetic  energy  « u^  + v*  + w* 

R Body  radius 

r Radial  coordinate;  density  ratio  o/p^ 

R Longitudinal  radius  of  curvature 

RF  Recovery  factor 

c 5 

R.  Reynolds  number  R.  » -2— 

^ 5 Vj 

T Temperature  (absolute) 

u Velocity  parallel  to  body  surface 

u Friction  velocity  u » / t /p 

X T W W 

Ug  "Wake  velocity" 

V Velocity  perpendicular  to  body  surface 

V Velocity  component  in  radial  direction 
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Velocity  component  parallel  to  axis  of  symmetry 

V Total  velocity  V = + v^  = ^ 

X Coordinate  parallel  to  body  surface 

X. j  Influence  coefficient  matrix  for 

y Coordinate  perpendicular  to  body  surface 

y*  y 

Y.  Influence  coefficient  matrix  for 

1 r 

z Coordinate  parallel  to  axis  of  symmetry 

a Body  slope  (angle) 

3 Boundary  layer  equilibrium  parameter 

V Specific  heat  ratio 

5 Boundary  layer  thickness 

5 Displacement  thickness  5*  s ( 1 - j ■ ) dy 
★ s s 

"Kinematic  displacement  thickness*' 

6*  = /^  (1  - u/Ug)dy 

£ Kinematic  eddy  viscosity 

Z Dummy  variable  of  integration 

n Isoparametric  element  coordinate 

8 Angle  of  velocity  vector  with  respect  to  z axis; 

momentum  thickness  9 = pu/p.u«  (1-u/u  )dy 

0 e e e 

0 Angle  of  velocity  vector  with  respect  to  body  surface 

k von  Karman  constant 

1/2 

X Skin  friction  parameter  X » C^/lC^| 

u Dynamic  viscosity 

V Kinematic  viscosity 

C y/d;  isoparametric  element  coordinate 


Dens  1 ty 

Surface  source  density 
Shear  stress 
Laminar  shear  stress 
Reynolds  shear  stress 

Velocity  ratio  u/u^ 

Turbulence  dissipation 

Finite  element  Interpolation  function 

Subscripts 

Edge  of  boundary  layer 
Reference  upstream  value 
eva 1 ua ted  at  wall 
Fluctuating  turbulent  quantity 
Stagnation;  Incompressible  flow 
■•Tall"  of  body 
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